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ABSTRACT 


Recent interest has centered on nonlinear techniques 
for the restoration of images. These techniques improve the 
resolution of blurred images without introducing erroneous 
details. Maximum Entropy method formulations incorporate 
these features and has been described and simulated on DEC- 1090 
computer system. Tlio iterative nature of the solution requires 
extensive use of the computer system. Restoration results 
are obtained for artificially blurred arbitrary image data. 

Minimum Entropy Deconvolution for restoration of images 
is another technique which is described and simulated for a 
class of signals which have essentially a non-Gaussian proba- 
bility distribution and are impulsive in nature. 

The description of various techniques and the software 

developed for their impler.ienta tion arc described in elaborate 
details. 
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CHAPTER 1 


INTRODUCTION 

The basic purpose of collecting data is to extract 
meaningful information about the phenomenon of interest. 

Often, the phenomenon is not a direct physical observable 
and the information about it is only ob'tainable through 
Some measuring instrument or system. The data available 
may be a linear superposition of the desired quantities. This 
type of linear information mixing is peculiar 'co a large 
number of diversified fields in physical sciences. The 
main problem confronted in these fields is hov/ to unmix the 
data, Viihen the phenomenon of interest is an image, the 
process of unmixing the data is termed as ’Image Restoration’ , 

Image restoration deals with images that have been 
recorded in the presence of one or more sources of degradation 
There are many sources of degradation in imaging systems. 

Some types of degradation affect only the gray levels of the 
individual image points, without introducing spatial blur. 

Such degradations are called point degradations. Other types, 
which involve blur, are called spatial degradations. 

For the purposes of digital Simula tion »based on the 

t 

convolutional nt>del , only spatial blur degradations win 
be considered. 
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In optical image formation, an unknov/n spatial radiance 
distribution o(x) called 'The Object' produces an image 
irradiance distribution i(y) which, for the purpose of digital 
processing, is collected as data i(ypj^)> (or for 

simplicity of notation as i(m), nt=l,.,.M), It is desired 
to estimate o(x) at a discrete subdivision of points , 
n = The relation between the image and the object, 

in one dimensional notation, is the linear form 
X 

i(Ym) = / o(x)h(y^,x)dx 4- vly^^), ntel,... k (l.l) 

where h(y,x) designates the point spread function of the 
overall image forming system and v denotes additive noise, 
random or otherwise. The assumption that noise is additive 
is subject to criti'^^m. However, the assumption of the 
additivity of noise rrekes the problem mathematically tractable 
The finite limit X denotes the finite extent of the iimaging 
system, i.e, the image forming system covers only a finite 
region from -X to X, such that 

h(y,x) = 0 for all [xj > X (1.2) 

The main interest is to estimate or restore o(x) only over 
this region. It has boon firmly established that the ability 
to restore an object image o(x),from the observation of its 
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degraded image i(yfjj)» is ultimately limited by the noise 
v(yjjj) in the acquired image data. More specifically, the 
ratio of signal to noise matters, and the mean squared 
restoration error is found to vary in an inverse 'nteiihcar 
with the signal to noise ration 

To facilitate finding o(x) at the discrete sub- 
divisionsx n = 1,...,N, it is convenient -to replace the 
integral in (l-l) by an approximating sum, by numerical 
quarature, such that 

N 

i(yjjj) = £ 0)^ ®(^n^ ^ 1,...,M 

n=l 

( 1 . 3 ) 

Numbers are input Vi/eights for facilitating accuracy in 

the approximating sum. If, except for translation, the 
degraded image of a point is independent of the position of 
the point, then the point spread function takes the form 
h(y-x) and eq, (l,l) becomes 
X 

= •/’ o(x)h(yjn-x^)dx + vCy^^), m=l,,..M (1.4) 

~X 

In this case the degradation is termed shift invariant. The 
assumption of shift invariance shall be considered throughout 
this thesis. 

From mathematical stand point, given the model (1.4) 
and the degraded image i(Yjyj)> the aim of restoration is to 
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make as good an estimate as possible of the original image 
o(x^). Evidently, any such estimation procedure must 
require some form of knowledge concerning the degradation 
function h(y,x) in (1,4). This a-priori knowledge can be 
effectively used to place constraints on the solution of the 
restoration algorithm. One such constraint is non-negativity 
of the restored pixels of the object image. 

There is no single optimum method for estimating a 
general o(x). There is a best restoring method relative to 
the type of a-priori inforrration regarding object and the 
noise one might have at hand. Such types of prior information 
about object shape or object and noise statistic are 
important in fostering a choice of tho restoration method. 

Other factors, such as the amount of data to be processed and 
the computing cost also influence the choice, 

/ 

The various restoration techniques can be devided 
under two main categories 

a ) L in ea r me tho ds 

b) Non-linear methods 

The linear methods consist mainly of techniques of 
inverting the discrete imaging equation (1.4) or linearly 
inverting or filtering the continuous version (l.l). These 
methods can be analysed by analytical means and are 
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widely used due to their coriputa tional speed. In contrasty 
the non-iinear methods are difficult to analyse and are heavy 
in computational cost. Very often their implementation 
requires iterative schemes of computation. Thus, the basic 
question that arises is why use the non-linear methods. The 
fact underlining the efficiency of non-linear methods is 
that for the optical signal images the linear methods are not 
optimum because these images in general do not follow normal 
sta tis tics , being always positive quantities. Also, the 
resolution of a blurred image of an object can be improved 
by various linear techniques, but oscillations or ringing 
around sharp changes of intensity in the image often induce 
erroneous details at the same time. What is required is a safe 
inversion technique, giving restorations free from ringing 
so that confidence can be placed in the reality of all the 
extra details revealed by the restoration method. Recent 
interest has therefore centred on non-linear techniques 
which incorporate constraints to reduce the artifacts generated 
in the restoration. 

Maximum entropy restoration falls under the category 
of non-linear methods. It is a fundamental method for the 
solution of the inversion problem in imago restoration. 

Stated most briefly, it selects the probability distribution 
that has maximum entropy permitted by the information we have 
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The reason for the use of maximum entropy criterion for 
restoration is as follows. In constrained deconvolution 
first a roughness measure is constructed which is then 
minimized subject to some constraints. Clearly many other 
roughness measures are possible, whose minimization or 
maximization would cause a smoothing influence on the 
restored images. For example, consider the function 

H=:-2xlnx , (1.5) 

where x In x win always exist since the image gray levels 
are essentially non-negative numbers. Structure of (1.5) 
is similar to that of conventional entropy measure which is 
defined for a discrete random process and gives its 
information rate* Here the entropy is defined for a single 
deterministic non-negative function and serves as its 
roughness measure. Maximization of this entropy measure causes 
a smoothing influence on the restoration. This smoothing 
influence is both on the noise as well as the object image. 

Thus a trade-off is achieved between the contrast of the 
object image and the effect of the noise. In certain algo- 
rithms noise can be smoothened out more than the object 
imago, thus producing bettor restorations. Also, due to the 
logarithmic nature of the entropy measure, its solution is 
deemed be positive and hence the positivity constraint 
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is taken care-off inherently. This aspect scores over the 
linear methods of the restoration. Furthermore, the maximum 
entropy solution results in the least biased estimates and 
protects the restoration against spurious details for which 
there is no evidence in the observed data. 

Chapter 2 of this thesis discusses one such nBxi— 
mum entropy technique for the restoration of images. This 
is based on the algorithm developed by B.Fl. Frieden and 

is essentially a (M+l] dimensional search problem, where 
M is the number of pixels being restored. Due to the 
large dimension of the problem as well as the iterative 
nature of the scheme for obtaining the solution, the compu*^- 
tational time requirements are very largo for images 
depicted by larger arrays. 

Chapter 3 discusses another maximum entropy 
technique wherein the search is directed along one dimension 
only for obtaining the desired solution for the restoration. 
Though this method is also iterative, but being a one dimen- 
sional search problem, it requires much lessor computa- 
tional time for larger arrays as compared to (M+l) dimen- 
sional search algorithm. 

In contrast to the principle of maximum entropy 
restoration is the concept of 'minimum entropy deconvolution* 



8 


This method is widely used in the field of seismology, 
where for a given time series y, which is a filtered version 
of Vi/hite noise x such that 

y = f X , (l»6) 

the deconvolution problem is to find a filter b which 
deconvolves or recovers the series x from the observed 
series y, i.e, 

X = b -»• y (1^7) 

This technique can be adapted for the purpose of 
restoration of a class of image signals v\fhich are impulsive 
and essentially have a non-Gaussian probability distribu- 
tion. The underlying principle, on which this method is 
based, is that the minimization of the entropy of the 
estimated image data restores the contrast or order in 
the estimated image. This is so because a data with 
non-Gaussian probability distribution has much order which 
can be restored by ensuring the solution to be as much 
away from the Gaussian distribution, as much permitted 
by the observed data. Thus, we require defining a norm that 
measures the order in the desired signal data. 

Another feature of this method is that the a-priori infor- 
mation about the point spread function of tho degrading 
imaging system is not needed at all. Thus, this technique 
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can be literally termed as ’Minimum Entropy Blind Decon~ 
volution’. Chapter 4 discusses all about this technique. 

The techniques discussed in Chapters 2,3 and 4 have been 
simulated on the computer DEC 1090 system. The computer 
implementation and the results obtaiiaed thereof are discussed 
in Chapter 5, This is folloV'/ed by the concluding remarks in 
Chapter 6, Tlie characterization of the digital image and the 
finite area superposition operator is discussed at Appendix A. 
Tlae computer program listings as well as the simulation results 
obtained are attached as Appendix B, 



CHAPTER 2 


MAXIMUM EHTRDPY RES'IDRATION : ERIE DEN » S APPROACH 

2,1 MXIMUM ENTROPY PRINCIPLE; 

&itxopy of a continuous probability , 

function p(x) is defined as 

H = - / p(x)in p(x) dx (2.1) 

— .eo 

and for the discrete case as, 

N 

H = - 2 p. log p. (2.2) 

i=l ^ ^ 

The rationale for using the principle of maximum 
entropy is as follows, A problem similar to the inversion 
of the basic imaging eq. (l,l) is the inversion of 

CO 

q(m) = / dx s(m) p(x), nt=l,2,...M (2.3) 

•^CO 

for the unknown probability density p(x) given M of its 
moments q(m) , which corresponds to the input image data 
i(yjjj) t p(x) corresponds to o(x) and must be positive to 
represent a real probability density, and s(m) is some 
weighting factors. For this probability-estima tion problem, 
it has been shown that the 'least-biased* estimate p(x) is 
the unique function which has a maximum entropy 
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H = - / dx p(x) In p(x) (2.4) 

-»oo 

Subject to the data cqn. (2.3). Least biased means as 
equipxobable as the input moments allow. Thus, the 

maximum entropy condition exerts a smoothing influence on 
the estimates, and ensures an all—posi tive output. 

2.2 fviATHEMATlGAL FORMULATION 

Measures of entropy used in imago restoration axe 
many, and have been widelyused in different imago processing 
schemes. Friedon's entropy measure [l,2] is discussed hero. 

The statistical model for the object is assumed to 
be composed of discrete, n© thema tical grains of small 
intensity &(o) which are distributed over the object scone. 
The scene is sub— divided into cells centered on a subdivision 
of points , and the unknown object is assumed to have 

0 ^ grains in cell n. Thus, recalling the one dimensional 
imaging eqn. (1.2), 

o(x^) = AqJ n=l,...N (2.5) 

Let p^ represent the probability of a grain locating 
in cell n. If a largo number of grains arc distributed over 
the object, then by the law of large numbers 


^n " ^n^^T' 
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where 0^ is the total number of grains. in the object* 0^ 
is assumed to be known from the image data by the law of 
conservation of energy. 

The entropy, then, is 

N 

- Z Pn Pn = ^ (0j.^/0^)ln(0^/0^) 

n=:l n 

Substituting 0^ from eq, (2.5), it follows 

H = ~ £ (o(x^)/£^^ 0^)ln (0^/A^ 0^) (2.6) 

n 

By the definition of 0^, V'/e finally get, 

H = -a Z o(x )ln o(x ) -}• p 
n 

where a,p are constants. 

The principle of maximum entropy then becomes 

H = - Z o(x^)ln o(Xj^) = maximum (2,7) 

n 

where ^(x^^) i^ the estimate of o(x^). 

The Frieden* s entropy measures, developed above, is 
essentially a concave function. Thus it has a unique 
maxima and hence a unique solution. The absolute maximum 
H in fact results from the perfectly smooth estimate 



13 


b(x) = Constant, The input image data, however, acts as 
constraints which force fluctuations in o ( x) and hence 
departure from the absolute maximum H situation, 

2.3 FRIEDEN‘S MAXIMUM ENTROPY SOLUTION; 

In this approach, a weighted sum of the image entropy 
and the observation noise entropy is maximized. An unbiased 
or maximum entropy estimate of noise is also desired now. 
However, the noise in the image data can be both positive 
and negative, and the logarithm of a negative quantity is 
undefined. This difficulty can be overcome by applying 
additional a-priori knowledge into the problem and assuming 
that a constant B is known such that 

B = - most negative vCy^^) (2,8) 

where vCy^^) is the observation noise at point m,in eqj(l,l). 
Thus, we define a new set of noise values 

V(m) = v(m) + B, B ^ 0 (2.9) 

B is a positive constant of sufficient size that it cancels 
the most negative going values v(m) and thereby makes all 
V(m) > 0. But Precise knowledge of this number B is not 
available. However, if a-priori information about cr, the 
standard deviation of noise, is at hand then we can set 
B = -2a. 
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The quality of the restorations do not critically 
depend on the use of a value B precisely satisfying eq, 
(2.8), however the solutions are best v^hen eq. (2,8) is 
indeed satisfied. 


Since the object and the noise are independent arrays 
of numbers, the total entropy from both is the sum of the 
two entropies. The input data which constrains the 
maximum entropy solutions away from perfectly flat estimates 
is the image data and its total energy The is 

equated to the total object flux by conservation of energy. 
In the digital representation of the image, implies the 
sum of the gray levels. The restoring algorithm, for one 
dimensional images, then becomes 

IJ ^ M A 

-2 o(x^) in 'o(x„)-f 2 V(y^)ln V(yJ 

n=l nfc=l 

- 2 /\(tii) [ 2 o(x_^)h(y_^ - + v(yJ-B-i(Yj] 

m=l n=l 

\ ^ A 

-■^(Mfl) ( E ~ I ) = maximam (2.10) 

n=l 


The new input parameter^ allows for affecting 
varying degree of smoothness on |^o(x^)j-and 

larger is the value of J^, the more emphasis is placed upon 
maximizing the entropy of the hence upon 
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making it smoother* The unknown ^ axe 

the Lagranges multipliers. 

For the two dimensional images the single summations 
are replaced by double summations and the single indices by 
double indices i,e, 

N N M M ^ A 

Z Z o(m,n) In o(ra,n) -EE V(m,n)ln V(m,n) 

nfc=l n=l m=l n=l 

MM. N N 

-E E A(p,q) S S o(m,n)h(p~m,q~n)+V(m,n)-B--i(p,q) 
p=l q=l in=i n=l 

. 0 N N A 

-/A(M'‘^+1) ( E 2 o(ra,n) - I ) = maximum (2.11) 

m=l n»l ® 

If the object and the image data are represented in 

vector-space form, then the two dimensional algorithm gets 

converted into the form (2.10) v^/ith numbers N and M 

9 9 

replaced by numbers N' and yr respectively, assuming that 
object and image data matrices were square. Thus, for 
analytical purposes, the objective function (2.10) will 
only be considered. 

The solution to (2.10) is obtained by differentia- 

/A ^ 

ting it with respect to o(j^j^) and V(y^) over all m,n and 
equating the result to zero. Differentiating (2.10) with 
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respect to o(x^) yields 

. M . 

-In o(x )-l - I /^(m) hCy^jj-x^)- /.(K'S+l) = 0 
nfcrl 

ox, 

M 

o(x ) = exp [-1- >.(1^-1) - E A (m)h(yj^-x^)] (2.12) 

ift=l 

A 

Differentiating (2,10) w.r.t. V(m) , we get 


-/in V(y^) -/~A(m) = 0 
or 

V(ym) = exp (-1- A(m)// ) (2.l3) 

The unknown lagrange multipliers (X( m)^ in the 
above solution are found by making the solutions (2.12) 
and (2,13) consistent with the input constraint equations 


S o(x )h(yjjpXj^)+V(yj^)-B = x(yj^); m=l,2,...M (2,14) 

n=l 

and, 

o(Xn) = Iq (2.15) 


The right-hand sides of (2.14) and (2,15) are the 
input data. Substituting (2.12) and (2.13) in (2.14) yields 

N M , 

E exp[-l- A(IVS+l)- £ /\('<)h(y.,x^)]h(y^,x^) + 
n=l ^=1 


exp[-l- A(m)/f]-B = i(Yjjj) 


( 2 . 16 ) 
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Eq« (2.16) represents a system of M non-linear 
equations in M+1 unknown , Mong vdth eqn. (2.15), 

it forms a system of Uhl equations, that are non-linear in 
M+1 unknown|A(in)J- , This system of equations can be solved 
to produce the desired restoration. 

There are many ways of finding the solution to the 
set of non-linear equations represented by (2.15) and (2.16), 
However, Newton-Raphson iterative scheme has been found to 
work effectively. In this, an initial set chosen 

for m = 1,2,. ..M and A(1 \/h- 1) is chosen so as to satisfy (j\h-l)th 
constraint equation (2,15). Usually X(m) is chosen = 0 for 
m=l,,,.M, The {"^mj are then changed in order to satisfy all 
the constraint equations. After 8 to -C) iterations of this 
algorithm, a set of rX(m)^ is found satisfying (2,15) and 
( 2 . 16 ). 


It can be seen that the functional form of represen ta 
tion (2.12) of the object guarantees a positive estimate. 
Regarding spurious oscillations in the output 0(x), we sec 
from (2.12) that 

d*^ o(x)/dx^ a o(x), k =1,2,... (2.17) 

This implies that o(x) cannot oscillate where it is 
at zero background and can oscillate little where o(x) is 
small. Hence, it is very smooth over these regions and 
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completely lacking in spurious details. Due to this very 
reason, this method works quite well on objects consisting 
of multiple impulses, for here all but a finite number of 
object values are zero. The physical objects consisting of 
rectangular pulses of random spacing and intensity are line 
spectra and star fields. 

The partial derivatives of o(x) with respect to X > 
m = 1,...,M+1, are also proportional to o(x), i.e. 




1 * • • ^ M+1 


a o(x) 


( 2 . 18 ) 


Therefore, vihexe o(x) approximates zero it v\/ill be 
nearly invariant with changes to the solution setJA^ . 

Hence, if the same object is restored many times by use of 
differing sets of noise-prone data, the solution 

I for these restorations v»fiii all be nearly invariant 
where o(x) v-r o. This amounts to a condition of stability 
to inputs where o(x)id^o, and hence is a very useful property. 


Further, the examination of object restoration eq.(2,12) 
reveals that each image input contributes a degree of 

freedom to o(x^). As the number of input image data decreases, 
o(x) approaches an increasingly smooth estimate. Conversely, 
larger M, freer is o(x) to fluctuate. 



CHAPTER 3 


MAXIMUM H’^TROPY EESTD RATION »ONE DIMENSIONAL SEARCH* 

3.1 INTRODUCTION: 

This chapter deals with an algorithm which involves 
solving a system of ordinary differential equations with 
appropriate initial conditions for the solution of maximum 
entropy restoration problem. This algorithm, developed by 
Zhung Yu and Haralick [17], ctoes not make use of any con- 
ventional optimization method. The main feature of this scheme 
is, that, it avoids searching in a (N+l) dimensional space as 
required by most maximum entropy algorithms, where N represents 
the number of pixels to be reconstructed in the desired 
image. This algorithm involves searching along a path in the 
(N+l) dimensional space which is defined by a differential 
equation system with appropriate initial condition which can 
be computed easily. Hence, this algorithm is essentially a 
one— dimensional search problem as compared to (N+l)— dimensional 
search problem in most ME algorithms, 

3.2 MATHE1\AATICAL FORMULATION: 

Let the underlying image be denoted by a set of 
positive numbers Defining the entropy measure on 

this set of numbers we get, 
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N 

H(Pj^ ... P^) = - Pi (3.1) 


where P. = f./ Z f., and N is the total number of pixels 

1 1 i 1 ■ 

being reconstructed* The observed data 
given by the convolutional sum 
N 

d. = 2 A., f. + e., j = 1,2,. ..M (3.2) 

3 i=l 31 1 3 


where e. *s are the additive noise terms with distribution 

given by N(0,a-), and A., represents the point spread function 
J J 

of the imaging system, Aji is a matrix of size MxN and is 
the convolution matrix of the imaging eq, (3,2), 

Let us define the following quadratic form which is 
the sum of weighted residuals 


Q( y ) 


1 ^ 
3=1 


1 ^ 


^3 


i=l 




A 


dj)' 


(3.3) 


The equation (3,3) can be expressed in the vector 
notation as follows 


Q(f ) = k (M. ■ diag 

Where, 

A = l^^ji^MxN 
£ = 


(3.4) 
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and D =; diag > . - . , 3 . 

For large values of M, the equation (3,3) can be 
replaced with the single constraint 

Q( ^1 » • « * » 

by Lavvf of Large Numbers, The condition Q = M determines 
the set of feasible images f, which passes the given 
statistical test for consistency with the measured data 
(3.2). 

Among the set of feasible images, the maximum entropy 
criterion selects the one which has the least configurational 
informa tion. 

Thus, we maximize the entropy subject 

to the constraints Q(f^, ,, , ,fj^) = “ M and Zfj^ = T, a const- 
ant, The Second constraint is introduced since the total 
intensity has a status different from individual pixel values 
and ensure the law of conservation of energy. Rewriting (3,1) 
we get 

N f f . 

H(Pj_,.,.,.Pj^) = - 3-og 277 

N f- 

= - 2 (log f^ - log Ef^) 

1— i 3, 
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:• N f. N f. 

= - Z log f. + log If. 2 

. = H(f^,...,fj^)/Ef^ + log 2fj^ 

From above , since Zf^ is a constant, it can be 
proved easily that the following three formulations of the 
. entropy measure aie equivalent to each other, 

a) Maximize H(Pj. , . . * ,Pj^) w,r,t, ( f , , . , , f^^j) , subject to 

constraints Q = ^ M and = T (3,6a) 

b) Maximize * . ,fj^) = -Ef^ log f^ w.r. t, f^, i=l. ..jN, 

under the constraints Q= M (3,6b) 

and = T 

c) Maximize [H( f^, , . . f^p+ix Zf^ - AQ(f^^ . . . ,fj^)] subject to 

Q = “ M and Z f . = T , (3.6c) 

^ i ^ 

where y, and /j are the Lagrange multipliers, 

\ie shall solve for case (3,6c) but without the 
constraint Zf^^ = T i.e. 

Max. [H(fj_, . . .fj^)+lJiZf^~ /^Q(fj_, . ,fj,^] subject to Q= ^ M. 


(3,7) 
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The parameter \x can always be adjusted such that the 
solution of (3.7) satisfies the constraint of total 
intensity. 

If f*, ..v, fj^ is a Solution of (3,7), then it also 

a solution of problems (3.6a), (3.6b) and (3.6c) when 

S f. = T. 
i ^ 

This algorithm for maximization of the entropy measure 
(3.7) is developed by first considering the unconstrained 
maximization problem, related to (3,7), as follows 


Maximize £H( fj^, , . , fj^)+ p.Zfj|^-" Q(fj^> •* »fjsj)l 

(= J(£; )) (3.8) 

Differentiating (3.8) with respect to f gives 


VJ 


dJ ^ ^ 
df df 


N 

[- 2 f^ log f. + f. -/Q(f, 

i=l ^ ^ i ^ 



= - [log fj. 




From (3.4) 

Q = 5 (M - d)^ 0 (^ - ^ 

= I A f - d’^D A f - fWs d + d”^ D dj 


(3.9) 
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^=|[2A'^DAf - A A - D d 3 

= I *[2 A^D(M ~ d )] ^''A’^D d = DA I 

= A^D (Af - d) 

Substituting this in (3,9a) vve get 

J = -[log %]+ - A,aJd(M-<^)- 

Equating VJ = 0 yields the solution of (3.8) 

Thus, if >0 is a solution of (3.8) then 

VJ = -[log f5,..,,log f^ ]+ (vi-l)[l,...,l]^^ 

A (Aff - d) = 0 (3.9) 

When /\2.0» the Jacobian of J 

y^J = -diag ^ ]- AA^D A < 0 (3.10) 

■^o — — - 

is negative definite. Hence J is a strictly a' concave function. 
Lemma 3 ,1. 

Assuming A >. 0, the solution of (3,9) could have at 
most one solution, and ’the solution (if it exists) must be tht 
maximal point of the function J(f; |i,/\). 

Proof ; For a strictly concave function, there exists at mc^ t 
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one maximal point. Since < o, hence S/J = 0 is the 

sufficient condition fox a maximal point, Alsoj» it is clea;, 
that when problem (3,8) is equivalent to solving the 

stationary point equation VJ == 0. 

Theorem 3,1 

Assame that > 0 and f° is the solution of 
f rl-i-y Aq ) = 0* Then p, AqJ uniquelly determine a 
branch f( ; p, f®) of maximal points of the function 

) corresponding to various values of A > 0. This 
branch of maximal points is governed by a initial value 
problem of differential equation, given by 

) fj = 

Aq, f®) = f (3.11) 

In expanded form, it becomes 

(diag ^ S A) ^ = - a’^D (^ ~ d) 

K A^; A^, f) = f (3.12) 

Proof; For Aq = o, if f° is the solution then 

Aq) = 0 and <^^3 if; p, A^) < O 

This implies that the initial value problem (3.11) has unique 
•■'Solution curve f(A ,p, whereAQ,^ belong to an open. 
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interval [ A 2 ) where A >_ 0 and /i 2 ~ ^ong this 
solution curve the gradient of S7 J with respect to A » 


dA 


v' J = 0. Since >7 J is a function of A and f, hence 


a r-r-r U r-r t / \ d'^/J df , U v U 

Tjr = "W" • dA ^ 


dVJ 


2 , df d -r 

= V cC^ + 35: ^ 

Substituting the value of n 7 J in the above equation gives 


VJ = (V^J) 37v+ ^ [-(log fj^,...,log + 

(p- 1 ) [l,,...!]^ -aaJd (M - d)3 
= iV^J) If - A^DCAf-d) 

or, 

^ '\/J = If . - \7Q(f) = 0 

Integrating from A^ to A » ''''® 9®"^ 

AA f P-f Aq » ^ ) f P'^A) 

u » A q) = 0 

l.e. VJ(f( >v ; ^ = V Aq) 

since i(^ 

But, f° being the solution of (3,8) for X - /<> . \7 J(f^ )-o 

. f H, f°); n,A) = 0 
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As a result, ; |a, a^, f°) is a branch of solution 

of the stationary point problem (3,8) or a branch of maxima} 
points of the function J(f; ) by Lemma 3,1, 

Lemma 3,2: 


Assume that a branch f(/\ } p,,,A^,f^) is defined by 
(3,12) where y J(f^ ^ = 0 as before. If there exists 

a point f = £( ^ ; Pj ^iong the branch such that 

Q(f*) = 0 (3.13) 

then the whole branch shrinks to a single point f® and for 
each i*. f^ = exp(ii~i). 

Proof ; 

Along the branch f(A ; \x, )\^,f°), we have 


^7 J(f(A ; n, n,A ) = 0 

OX 

-[log f^,.,.,iog 3^ + (p-i)[l,.,,,l]'^ =XVQ(f) 


From this it is clear that 

Q(f ) =0 implies f « ».« 
(3.11) is equivalent to 


^ = VQ<f) 



* 

A= A 


* 

f 


fj^ = exp (p-l) and 
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The initial value problem thus has a unique constant 

* 

solution given by f = f * 

Theorem 3,2 ? 

Assume that ^0 and ^ is the solution of 
J(f; = 0- If for each S, 




^IS 

1 

r‘ - 

0 

• 

» 


• ; 

• 

! 

^US 


• 

0 





f S = 1 y * * *N 


(3.14) 


then the branch f(.A ; p-# /. q» f^) could be extended to [o,«). 


Proof; 


Let the branch be defined on A 2 ^ where A 0 


and A 2 ^ Along the branch we have 
V J = 0 ' 

i • ^ • 

"•(lo g f • ^log fjijjJ +’()j'”"l) ~ AiiSA^. 

- A d (3,14a) 

Le-t P be A D A . From assumption (3.14), it is easy to 
see that 

.2 


Pcc = 2 

K=1 o, 


^ 0 } S = 1, .,.,N 


K 


V-ie could prove that f(A ; p, A^jf^) is bounded from 
above when A~^A 2 » Assume the contrary,, i.e. if there is a 
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sequence A that 

(A^; f^) = max f°)--^ 


as K — ^ CO 


then from (3,14a), by dividing each term by f , we have 


log f 


^ , (u-1) _ 


f . 








^ +A,. P, 




The left hand side tends to zero and the right hand 


side would be greater than some positive number because 
Pgg > 0 and A 2 ^ ® v\tien is closed to A 2 * This being ? 
condradic tion, we conclude that there exists a positive 
number L > 0 such that for all A £ 

fj^CA ; p, A^f f®) < L ; i = 1,...N (3.14b) 

If A 2 ^ '"y then (3,14b) would guarantee that the 
left hand side of (3,14a) is bounded since the right hand 
side is bounded whenAE [Ajl» ^ 2 ^* this case, there also 
exists a positive number g> 0 such that VAE[A 3 ^» 

A ^ fJ'f Aq> ^ ^ B • (3,l4c) 

Because of (3.i4b) and (3,l4c) for any sequence Aj<^ 
there exists a subsequence and such that 
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Therefore 

^7 . . ,fj^; lir A 2 ) =0 and f'^ is the unique maximal 

point of J(f^p,A2) Lemma 3,1. Since any convergent 
subsequence f( Ak » Aj.? tends to 

.JT) 

the same limit point f 3 the sequence ^(A k' 
itself would also tend to f . Therefore 

lim f( A? ~ -* 

A-^A2 

Thus f ( A ? Pj /\q> can be extended beyond A This 

leads to the final limit A 2 = “• Similar arguments leads to 

~ 0 . 

Cox ; 

Under the assumptioiis of theorem 3,2, there exists a 
positive constant L such that [o,**), 

qC A; ii.Ao, t°) < L; 3u=i,...N 

Remark ; 

If for some S 



then, P» A ) =0 if and only if 

V i ^ 
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d . 

- log f. + (tJi-l)=: E P.. f. - S A.^ “x-l 

1 j^S ij 1 j ij Q-j 

Now fg is determined from ^s+l^*''*^* 

Corresponding to (3.15), A' = [A^, . , . . . ,Aj^] . 

Thus, the condition (3, 14) is not a troublesome restriction 

3.3 SOLUTION OF DIFFERENTIAL EQUATIONS: 

Vi/hen^^ = 0, the solution of the stationary point 
eqn. (3,12) is 

f° = f® = ... = :^ = exp (ii-l) (3.16) 

According to lanma 3.2, the following initial 
value problem 

y^J(f; p, ;\ ) ^ = ‘7 Q(f) (3.17) 

f] = [Gxp(p-l) , . . V ,exp(p-l)^ 

A= o 

determines a branch f(^ , p,/iO,i°) of maximal points of 
•J(f;p, X ) for 0 < X < «> When V s , Ag 0 

If it happens that Q(exp(p-l) , , . . ,exp(p-l) ) = ^ then 
f^= f^= ... = :^= exp(p-l) is exactly the solution of 
problem (3.6c), and it’s not necessary to solve the initial 
value problem (3.17). If it is not so, then we have to choose 
p in such a way so that a point f^,...,f^ on the branch 
f( A ; p,o,f°) would satisfy the condition Q(f ) = 2 M, 
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The solution curve of the initial value problem (3.17), 

corresponding to various values of \x covers the maximal point 

set of J(fjj^, . . , ,fj^; Vi»A) with respect to various p and A (A > o) . 

Moreover, they coincide with each other. Our problem is to 

choose ]x in order that the co.rresponding branch f(>, j p, o, 

nr 

[exp(p,~l),.,.exp(jA-l)] ) will now meet the condition 

q( f , • ■ • , fj^ ^ ~ 5" 


Theorem 3,3: 


Assuming that (3,14) holds, then for o < A < either 
t o, [exp(ii~l) , . ,exp(p~l)] ) = exp(p.~l); i=:l,...N 
or 

j T" 

^ Q(f(A ; ix,o,[exp(^~l),..,exp(ti-l)] )) < o 

(3.18) 

Proof ; 

The branch f(A ; m o,f°) given by initial value 
problem (3,17) is defined on [o,«j) by theorem 3,2. Along 
the branch we derive 


dX ~ ^ • 


df 

dX 


,T df 


‘A 


= (yQ) 

H ’f’ 

Substituting the value of ^ from (3,11), we get 

^ = (VQ)^ (V Q) 
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Now the following two cases arise; 

a) For some o < )\' < «>, y Q(f( A i.yojf'^)) = 0. 

In this case, fox o < }\ < 

r P-» o,f®) = exp(p,~l) by lemma 3,2, 

b) 'y Q o along the branch. In this case, < o since 

2 

J < o. Hence the proof. 

Let us define sets U,V,v7 as follows; 

U Q(exp(p-l),. . ,exp(p.-l)) = 

V Q(exp(p-l),.,.,exp(p-l))> M } 

V7 = f|if Q(exp(p-l),.,*. ,exp(p-l)) Or 

{3.19) 

It is amply clear from theorem 3,3 that only those 
initial value problems would produce meaningful result whose 
p belongs either 'bo set U or to set (vrWJ). In other words, 
only those values of p are useful which results in f ( /\ ; p,o,f*^) 
meeting the condition Q(f) = 

If pt U, then f^ = f 2 = % = expCp-l) is the 

solution to problem (3,7) with = o. It is also the 

N 

solution of problem (3.6a) and (3,6b) with E exp(p-^l) = T, 

i=l 

If p4 ivAwl , then Q(f (o ;p ,o, [exp(p-l) , , . , , exp(p-l)] ’^) 

1 ^ t J 

> M and Q(f( X ;p,o,[exp(p-l) , , . . , exp(p-l)] ) decreases 
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strictly nonotonously as _/\ increases from 0. Somewhere f 
Would reach a value such that Q(f) = ^ M, Thus, if for 

some A= A r Q(f( A "f t » • • • )= 2^^’ 

f( [exp(ti-l),. . .exp(p,-l)]^) maximizes the entropy 

H(P^, ,Pj^) subject 'to the constraints 

Q = h M and 2 f . = E f* . 

^ i i 

Computing sets U.V/W 

A direct computation leads to a quadratic form 


q( a,..*, a) = -h ba + c 


(3.20) 


The quadratic form (3,4) yields 

N 


a 1 I 


(3.20a) 


N 

M 

b = - Z — 

j=l a2 


(3.20b) 


d? 


c = §• I 

^3 


(3.20c) 


It is clear that 


b < 4ac 


Set U 


Vjhile computing the set U, the following cases may arise 
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i) If a > o, b^-4a(c - ^ lA) > Ot c>|'M and b < o, then 
Q(a,.,.a) - ^ M = o has two positive roots and a 2 » Thus , 
the two values of p, constituting the st U are 

exp 

or - 1 + log 

and similarly 

^2 = 1 + log a 2 

Hence US ^1+ log 1 + log (X 2 | 

ii) If a>o, c = |- M and - > o, then Q(oc,. . » ,a)- M = o 

has one positive root and one negative root. Since log 
of negative number is undefined, so we consider only the 
positive root. 

Thus, U= J 1 + log 

iii) If a > 0 , c < ^ M, it has one positive root and 
Us i 1 + log a, i , 

I v> 

iv) If a = o and c = ^ M, then U is an infinite open 
interval i,e. 

U = 00 )j , 

Apart from above four cases, U = 0, an empty set. Also 
the fourth case occurs very rarely. In general, U is empty 
or has two elements at most. This corresponds to no equi- 
entropy solution or two equi— entropy solutions. 
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Set V i 

V = ] p,; Q(exp(p-l) , , ,exp(p-l)) > h u\ 

i . ■"’ 

Here also, while computing the set V, the following 
cases may arise: 

i) If a > o and b - 4 a (c - 2 then due to the 

inequality of the quadratic eqn. Q(a,*..,a )-~M>o, V 
is an infinite open interval 

V » f ( — oo, «>) '• ( 3 . 20a ) 

J 

r\ 1 

ii) If a > o and b - 4 a(c - ^ M) > o, then Q(Q:,...a)- 
^ M = o has two real roots a^, and due to inequality, 

the set V will consist of two open intervals* 

Vr. 1 + log(max(o,a^))), (l+log(max(o ,0^) ) ,~)|- 

( 3 . 20 b) 

where ^ * *^2^ * 0(^2? * • # ,^2^ “ ^ 

iii) If a = o and b > o, then 

. ^ M-c > 

V =*i^l+log (max( 0, g— )),. “j 

iv) If a = b = 0, and when c > ^’ M, then V= (~oo,oo) - 
and vjhen c ^ ^ M, Y ■= 0 
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v) If a s» 0 and b < o, then 


V = ) 1 + log max (o, )) 


The case (iv) occurs rarely. In general,. V consists 
of atmost two open interval. 

Set W ; 

V/=|p.; ^Q(€xp(p.- 1 ) , . . . ,exp(.p-l)) ^ o’ 

Substituting the value of y Q? can write 




ip; exp(p-l)A^D A A^D d 9 ^ 

L ■ LiJ " ^ 

Vihile Computing W, the following cases may arise:- 



r. 

1 

i) If rank of X = 

ADA 

• 



^ T 


(3.20c) 


then W = )J 

ii) If rank of the matrix defined above = 1 and 


Q 6. • ^ Of then there exists a unique a such 

Li-j 


tha t 


a a"^ D A : d 


(3.20d) 



38 


iii) If rank of X is 1 and also 


A^D A 



0 , then W = 0 


In image restoration problems, the case (iii) occurs 
very rarely and generally the problem falls under category (i)* 

3.4 DEVELOPMe^T OF THE ALGORITHM: 


In practical image restoration problem vve never 
compute the set U as the solution so obtained, though it 
maximizes the entropy absolutely, is completely devoid 
of any contrasting features and depicts a uniformly smooth 
background. Hence, we have to look for a iji belonging to the 



Rewriting the initial value differential equation 

(3.12),: 


(diag ]+ A^D A ) ^ = -A^D (M “ ^ ) 

f(o; \i,f o, f°) = f° = [exp(p~l), ...,exp(^~l)]’^ 

From this differential equation systeni,v;e can easily 
obtain the following iterative equations as follows; 

(£ + a’^D a ) df = -dA - d ) 

or in the form of an iterative scheme. 



39 


(£^ + A^D A) - f^) 


-(Ak+l-Ak) a'^D(AX~ d), k = 0,1,2,... (3.21) 


where. 


F — diag [ , » * • , , ] 




^ ° ’ A k+i > 


V-?hen each 1 A k+i small enough, each would 

approxin^ite o»f°) very well. Rewriting (3.21) 

as follows 


(e’^ +Aic a'^S a) = (2Aij -A^+i) a'^d a 


+ [1...!]'’'+ A^DA (3.22) 

The expression (3.22) represents a large system of 
linear equations. It can be seen that the structure of 
the coefficient matrix (F^ k— such that it is 
positive definite. 

The above system of linear equations can be solved by 
various iterative schemes, but due to the structure of die 
coefficient matrix. Gauss- Seidal iterative scheme can be 
employed to solve (3.22) efficiently. The Gc. us s- Seidal method 
[21] is a basic scheme for solving a system of 
partial differential equations iteratively. 
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Let g_ represent the coefficient matrix (F^ k Q 
Let b represent the vector; 

b = (2>i}^ - A^D A f 4 - [1,..,,!]^ + 

+ (/N k+1 d 

Setting up the iterative equation for 


f^= exp (p-l) 


1 

« 

1 


:k+l 


^ii ^ 


j<i ^ 


k+1 


Z P. .f^^ ] 
j>i ^ 


for 




(3.23) 


Since 1 j\ “-^k^ appropriately small, f^”^^ shoul 


d 


be near to f^. This greatly speeds up the convergence of (3,23) 


To summarize the algorithm, a flow diagram is presented 
as shown in Fig, 3.1 


From the flow chart it is clear that, having fixed 
the value p, we are only varying the value of to achieve the 
desired solution. Hence this algorithm is essentially an one 
dimensional search problem. Due to this reason, the algorithm 
is more efficient in terms of computing time as compared to 
the earlier (N+l) dimensional search algorithm. However, the 
restorations obtained by this algorithm are somevjha t inferior. 
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Another algorithm which is faster and more efficient 
than the (N+l) dimensional search algorithm and claims to 
produce better restorations than the one-dimensional search 
algorithm ,has been developed by BURCH, GULL and SKILLING[19], 

A brief outline of this pov/erful maximum entropy method is 
discussed below, though its implementation is beyond the 
perview of this thesis, 

3.5 POV’/ERFUL MXIMUM ENTROPY METHOD: 

This method is applicable to a wide class of problems 
including that of image restoration vihen the data are the 
convolution of an original object v/ith a spatially invariant 
point-spread function, with the addition of random noise. 

The required restoration is a set of positive numbers 

f ,f, ,,.,f., T which are to be determined, and on which the 
o 1 N— 1 

entropy 

N-1 

S = - Z p'(j) log p(j) p(j)= f(j)/i:f (3.24) 

j=o 

is defined. Related to the set of numbers f(j) is another 
set of numbers g(i) vjhich represent the simulated data and which 
would be produced (in the absence of noise) by the observing 
equipment if the original object were correctly represented 
by the numbers f. Generally there can be an arbitrary 
relation between g and f, but the here we will consider the 
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case of convolution with a known spatially invariant p.s.f. 
or blurring function h. Thus 

N-1 

g(i) = X f(j) h(i-j) ; i = 0,1, ...N-1 
j=o 

where the suffices are written as one dimensional for 
simplicity. Setting up some statistic to measure the misfit 
between the actual noisy data d(i) and the simulated data g(i), 
the most convenient choice is Chi- Squared, 

(g(i) - d{±)f/ a^(i) (3.25) 

i=l 

2 

For large N, the plausible values of X are all close to 
2 2 2 

X„ = N. The condition X = X now determines the set of feasible 
o o 

images f, which passes the given statistical test for consis- 
tency with the actual data. Among the fasible images, 
maximum entropy criterion selects that particular f which has 
least configurational information. 

The maximum entropy solution can be obtained by 
maximizing S(f) over X^(f) = . The required maximum of S 

r\ 

will be at an extremum of Q = S- /\ X for sui'table lagrange 
multiplier A. , From this f can be determined via 

log A - log f(j) = (Ax f) dx^/df(j) (3,26) 

3 

where 


A = exp (x P(i) log f(i)) 
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The quantity A is a weighted mean of thef(i)and can be 
treated as the default value to which f(j) tends if there are 
no data pertaining to cell j (dx /6f(j) = O), 

In image restoration, since the total intensity 2f 
has a status different from individual pixel values, hence Sf 
forms a second constraint with its own multiplier p. Thus 
Q gets modified to Q = S - P-Sfy which in turn modifies 

K to 

A = exp (Z p(i) log f(i)-- p2f(i)) (3.27) 

One can either choose p and hence A, to fit given value 
of Ef or one can simply use A itself as a predetermined 
parameter. Generally, the latter choice results in faster 
implementation of the algorithm. This is equivalent to 
maximizing a modified form of entropy 

S = -2f( j)log(f( j)/ eA) (3.28) 

q 2 

over X = X , for which the external solution is 
o 

log A - log f(j) = dxVdf(j) (3.29) 

The results are not sensitive to the precise value of A, 



45 


Normally, such a problem is solved with the aid of a 
lagrange multiplier /v vjhich transforms it ixi an unconstrained 

O 

maximization of Q(f) = S r-A. X . The simplest such algorithm 
is steepest ascent, in which one performs a line search 
along the direction V Q- This is however inefficient. The 
standard way of improving a steepest ascent algorithm is to 
use some variant of the conjugate gradient technique. Instead 
of using •CyQ itself as the search direction, one uses a 
combination of this with previous gradients choosen to 
give quick convergence if Q is exactly quadratic in f. 

However, Q in maximum entropy is highly non-cjuadra tic and 
conjugate gradients also proved inadequate. Thus, to make 
the algorithm more powerful, instead of searching along just 
one direction at a time, a subspace of three directions 
was used, namely 


(ei)i = f(i) ds/df(i) 
(e^)^ = f(i}e)x2/df(i) 

.2v2 


(®3h = Sflfr 

where, in the definition of Q, is given by 
= (2f(iXb3/df(ih^/2 f( i)(c) X^/df( i) 

Thus a simultaneous search along three directions is carried out. 

Convergence to the image solution is tested by checking 

o 

that the directions 'vX and s axe parallel to within a 
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o 9 

few degrees and that X = with less than about 1 /. errox. 

The number of iterations needed depends on the signal to 
noise ratio, roughly as its square root, and 15— 2D itera- 
tions are generally sufficient. Fortunately, this number 
does not depend on the no', of pixels N being reconstructed 

So that overall computing time has N log N dependence. 

Approximate computing time for images of various 
sizes are given in Table 1 for various computer systems. 

Table l: ivlE Computing Time 


Programmed Actual Time normalized 



size 

CPU time 

to 512x512 

IBM 3033 

256x256 

8 min 

32 

min 

PDP 11/60 

512x512 

20 hr 

20 

hr 

PDP 11/34 

512x512 

50 hr 

80 

hr 

POP 11/34 + AP120B 

512x512 

2 hr 

2 

hr 

NORD 50/10 

1024x1024 

20 hr 

8 

hr 

VAX-11/780 

512x512 

3 hr 

5 

hr 


CHAPTER 4 


MINIMUM HnITROPY DECONVOLUTION 

4.1 INTRODUCTION: 

Minimum entropy deconvolution (MED) was developed, by 
Donoho [ 22 ] 3nd Wiggins [23] fox the extraction of infor- 
mation about the earth’s reflectivity function in the field 
of Seismology , The input signal in such cases is always 
consisting of few spikes and the. output is a time series 
whi<;Lh is the convolution of the input signal with the 
relectivity coefficients of the earth’s crest. Thus, the 
only a-priori information available is the general nature of 
the input signal. The idea of MED can be adapted for the 
purposes of digital image restoration for a class of signals 
which are impulsive and having a non-Gauss ian probability 
density function, . The mathematical justification of the 
principle however remains essentially same, 

4.2 MINIMUvl ENTROPY PRINCIPLE; 

vcts/on 

Given a time series y, which is a filtered^of a 
white noise x 

y = f * X (4.1) 

where f is the convolution filter, the deconvolution problem 
is to find a filter b which recovers x from the observed 
series y, such that 
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X = b * y (4,2) 

Since the delay characteristic of the convolution 
filter f is not generally known, the second-order methods 
(spectrum, correlation) are unable to solve the deconvolution 
problem in general. For a typical filter of length p+1 , the 
number of filters with distinct delay characteristics but 
same second-order property are 2^. 

Minimum entropy deconvolution is a technique for 
deconvolution without making prior assumption about the delay 
characteristic of the filter f. In this approach we seek 
an operator which, when convolved with the observed signal, 
converts it into a signal with ’simple' appearance. Here , 
the term simple implies that the desired signal, recovered 
from a signal smoothed by convolution with some system 
function, is impulsive. Such an approach maximizes the 
order or equivalently minimizes the entropy of the signals. 
Hence the name Minimum Entropy Deconvolution, 

The idea of simple appearance or maximum order or 
minimum entropy are illustrated in Fig, 4,1. It shows a 
time series generated according to (4.1), where f is a mixed 
delay filter. The result y is then filtered using b'*’, b"",- 
and b filters, all with same second order properties, and 
correspond to maximum, minimum and mixed delay assumption 
about f respectively. 
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tu It t<S 
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It is apparent that the recovered signal x corres- 
ponding to the correct assumption about the delay chara- 
cteristics of the filter f has the simplest appearance. 

Thus, the eye, using the judgement of simplicity, can identify 
the correct solution to the deconvolution problem even 
though the correlation and the spectrum methods could not. 

Thus we are looking for a function of the data that reproduces 
roughly the judgements of simplicity that the eye could make. 
If O is such an objective function, then our aim is to obtain 
a formal- estimate fox the deconvolution problem i.e, 

b - Maximizer cCb'^^-y) (4,3) 

b 

-'N ^ /V 

where b is the filter which gives x = b y the simplest 
appearance among all series of the form x = b * y. 

4.3 OBJECTIVE MOTIONS; 

The objective functions, which work very well for 
the deconvolution problem discussed above, are 


i) 



(4.4) 


where is the finite sampled version of the estimated 
series x. The objective function 



4 \ 

02(xp = 


(4,4a) 
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is a Special case of (4,4)* 

ii) The objective function based on an information theoretic 
point of view is , 

Op (x^) = -H ) (4*5) 

where H represents the entropy function and is a scaling 

factor. 

All these objective functions have two properties in 

common, 

a) They are scale invariant i.e, 

0(ax) = 0(x) for a o. 

b) They depend on the sample only as an unstructured 

batch of numbers i.e.^tl^ey depend only upon the empirical 
distribution of the number (x(l) , . * . ,x(n) ) viewed as a 
simple random sample. 

Due to above properties, 'the MED type procedure can 
be successful only if there is a difference between 'the 
empirical distribution of X and that of x = b y, where 

^ 5^ b. 

Due to the probabilistic structure of the problem, a 
partial order is introduced which will help in giving fairly 
complete characterization of the objective functions which 
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are admissible fox use in rule (4,3). For the purposes 
of analysis, we shall consider the input signal as a time 
series only. The partial order and the assumptions on the 
data series which allows this method to work are as follows; 

4.4 PARTIAL ORDER •>: 

The deconvolved time series x can be written as 
X =: "b * y. Substituting for y, 

X = b * f * X 
= (b f) ^ X 

Consequently, each series is sampled from random 
variables which are linear combinations of independent and 
identically distributed random variables i.e. is sampled 
from ^t-u * '"^here g = b*f* V/e shall be discussing the 
distributional properties of the linear combinations and their 
implications on model (4,l)-(4,2). 

.Let £, be a set of random variables with finite 
variances which is closed, under linear combinations, such that, 

2 a^P^ + c is in P for all P i^-E. 

Two random variables p and Q will be regarded as equivalent, 
written P = Q, if for some constants c and a o , ap+0 has 
the same probability distribution as Q» = is an equivalence 
relation on p , 
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Definition of Parial Orde r : 

P •_> Q implies that fox appropriate constants 

with 2 9? < “i 

Q = Z a .P . 

v^here the P^ are independent copies of P. The relation -> is 
short for ^ not i , 

’>. a partial order on i> because of the following 
tv/o properties 

a) Transitivity;- If P Q and Q R, then 

P •> R 

b) Asymmetry:- Let P and Q have finite variances. If 
P Q and Q P then P i Q, 

proof ; Property (a) follows immediately from definition, 
since if R = S Q = Sa^^P^, then 

■"SR- "1*^3 '’iJ 

r , j 

To prove property (b), the following characterization 
is required: 

R is a Gaussian random variable if and only if R 
has finite variance and 


R = S a.R. 


( 4 . 5 ) 
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holds for some coefficients v^hich are non trivial, i.e, 

at least two > o, and square summable < o°)» 

Now, suppose Q = X P = Xb^Q^. Then 

P = E a.b^P. .. This is precisely a relation of the form 

j J <1- ^ 

(4,65. So, if either or nontrivial., then P is 

Gaussian random variable. Thus, the characterization implies 
that any linear combination of P^s, e.g, Q, is also Gaussian, 
in Vihich case P = Q, If the coefficients 
both trivial, the Q = P is self evident. 

The characterization results provide that if R is 
Gaussian, then R P for any P, since any linear combination 
of copies of H is Gaussian, and hence equivale.nt to R, On the 
other hand, if P has zero mean and finite variance, then 
from central limit theorem, 

EPj|_ converges in distribution to a Gaussian 

random variable R* Since p > S and 3^--> R in distribution, we 

n n 

can conclude that p R for every P in ^ . Thus, is a 
partial order with this additional stipulation. Hence the proof. 

Thus, for P in P and R Gaussian, 

P '> 2 a^P^ •> R 

This order is strict unless either; 
i) P is Gaussian. If so, then P = ^®i^i ^ 


(4.7) 

(4,7) 


i} 


are 
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ii) P is not Gaussian, but the linear combination is trivial 
(no two a^‘s non-zero). In this case P i 

From relation (4,7), we derive that the linear combinations of 
independent random variables are more nearly Gaussian 'than the 
individual components of the combination. Thus, v/e introduce 
a mnemonic ’P •>_ Q* Implying 'Q is more Gaussian than P’ . 

The equation (4.7) has a straightforward translation 
to a time series problem, and we can say that every filtered 
version of a white noise is more nearly Gaussian than the 
white noise itself. For example, a white noise p is naturally 
associated with the random variable P from which is 
sampled. A filtered version of P, h''''p, is associated with 
Zh^Pt_u* Then for filtered white noises p and r, p •>. r is 
defined to mean P b R fo3: associated random variables. Then, 
(4,7) becomes 

p > h ^ p r (4,3) 

where p is white noise sampled from copies of ^ in £ , and r 
is Gaussian time series. 

In the deconvolution problem, given by (4.1) and (4.2), 
the above characterization has immediate application. Since 
X = b * y, and substituting h = b f in (4,8) we got 

b * y > b * y ->_ z 


(4.9) 
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where z is a Gaussian time series. Thus, filter b is 
characterized by giving the least Gaussian output of any 
filtered version of y. This characterization of b by (4.9) 
makes the numerical maxindza tion in rule (4,3) possible. 

It will be seen later that if the objective function 0 
agrees with the partial order then the rule (4,3) is a 
consistent way of estimating b, 

4.5 UNIQUENESS OF THE MODEL. 

The model, represented by convolutional expressions 
(4,l) and (4,2), is essentially unique when the representation 
y = f * X is satisfied only by f and x or by scaled/shifted 
versions of these. 

L amma 4, 1 ; The model 

y = f X , where (4,10) 

where x is a white noise sampled from x in X and f is an 
invertible filter, is essentially unique if and only if X is 
a non-Gaus sian random variable. 

Proof ; 

Let there be tv-^o representation satisfying (4.10), 

i. e, 

y = f ^ X, 


y = f* * x’ 



57 


Novi? f X and f* ^ x* have the same power spectrum, 
and both f and f invertible. Thus we write 

x ’ = ( f ’ ■ f ) * X 

X = (f-^ ^'5- f I ) x‘ 

Both X •>, Sind x' x i,e, x' = x. 

Now, from relation (4,8) it follows that either X is 
Gaussian, or f~^ * f is trivial. Hence for the model 
to be unique, X should be non-Gaussian and the two represen- 
tation can differ only in scaling and time origin. 

To illustrate the non-uniqueness of the convolution 
model, when x comes from a Gaussian process, see Fig, 4,2, 

The deconvolutions are affected as in Fig, 4,1, Vie can see 
that the three results x"^, x”" and x are all uncorrelated 
and statistically undistinguishable. Also, the eye can not 
see any qualitative difference between the results. Hence, 
the judgement of ‘simple appearance’ is practically impossible 
if the input signals were Gaussian, It can be now said that 
eye makes the same sort of qualitative distinctions as 
dictated by the partial order *>,• 

4,6 CR^RACTERIZ/^.TION OF AmiSSIBLE OBJECTIVES: 

To find out under v/hat condi'tions, on objective 
function 0, will the estimation rule (4*3) will be consistent, 
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it is essential to characterize the admissible objective 
function. The basic idea is as folio v/s. 

By the ergodic theorem, the impirical distribution 
of converges as n increases to the distribution of the 
random variable X associated with x. Now, is a good 
estimate of x to the same extent that X resembles X 
rather than a linear combination Sa^X^, So the objective 
function 0 has to discriminate between the distribution of 
X and that of Sa^X^, Thus, 0 must do roughly the same 
thing as the partial order •>, 

Def init ion ; Objective function 0 agrees with •> in X if 

X •> y implies 0(x) > 0(y) (4,11) 

for every x,y which are filtered white noises sampled from 
random variables in X. 

. The restrictions on 0 is that it must be scale 
invariant since the scale of x and b are unknown and must 
be continuous function of the distribution of. its argument. 

Cons is tency ; 

The estimation rule (4*3) for a finite length data 
series can be written as 

b^ = maximizer 0 (b * y'^) (4,12) 

b ; length (b) = p+1 
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If the objective function 0 is scale invariant and 
a continuous functional, then sequence (b^), def ined above, 
is a consistent sequence of estimates of b for every non- 
Gaussian y^of the form (4. l)^ if and only if 0 agrees with '> 
on X, 

This shows that an objective which agrees vdth •> 
will give consistent estimates no matter what the distri- 
bution of X, as long as it is non-Gaussian. 

4.7 IMAGE RESTORATION PROBLEM: 

The model described so far can be fitted in many 
fields related to signal processing. In the field of image 
processing, the equation (4,1) describes the blurring/degra- 
dation of the object image and equation (4,2) describes the 
deblurring or deconvolution operation for the restoration of 
the object image. 

In this setting, a picture f which has been blurred 
by an unknown filter h can be recovered, provided the inten- 
sity or gray level of the pixels f does not follow a 
Gaussian distribution. 

Let f represent the digital image vector 
f = f ^ , , . , , fj^ 

and h represents the impulse response coefficient vector 
of the imaging system. 


!l — y • • • > 
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The degraded/ blur red image vector d is given by 
d. = I f. j = (4.12) 

j i ^ 

V'lhexe M = N+L~l 

The vector d represents the observation data* 

In the absence of the a-priori information about the 
p»s.f# of the imaging system, the problem of deconvolution 
or restoration is formidable* But this problem can be 
tackled when it is known that the object image has a gray 
level distribution v;hich is non~Gaussian and the gray level 
variations are impulse like. 

In the deconvolution problem w'e are looking for an 
inverse filter g of length K which, when convolved with the 
object image , minimizes the entropy of the estimated 

data. The order or length of the deconvolution filter has 
to be small as compared to the length of the observed data 
because the discrete convolution results in spreading of 
the data. The resulting deconvolution equation is thus 

^ . jM+K-1 

f = 2 d. . g. ; j = 1,...,M+K-1 (4.13) 

^ i=l 

where K is the order of the deconvolution filter, and f is 
the estimate of object image f. 
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4.8 SOLUTION; 

The Solution of the posed deconvolution problem rests 
on defining an objective function that measures the simpli- 
city of the desired signal. In other words, we seek a 
norm v/hich minimizes the entropy such that the order is 
restored in the estimated data. 

The two such objective functions, as discussed earlier 
and defined by expressions(4,4) and (4.5) are 

i) Oy. E /( E f (4.14) 

3 3 

where f^ is given by (4,13) 
ii) Og = “[- S fj log fj] = “H(fj) 

= S f • log f • (4,15) 

j J J 

where H(f .) is a measure of entropy in the estimate f. 

V/ 

The objective (i) is called Varimax norm, whereas 
(ii) is called an entropy measure, ^Maximization of those 
objective results in minimization of entropy i,e, maxi- 
mization of the order in the desired signal. The minimization 
problem is not. an unconstrained problem because the image 
gray levels are always positive and also the total intensity 
has a status which is different from the intensity 
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of individual pixel. Thus, the problem develops into 

Maximize (0^ ox Og) subject to the constraints 

£ 


a) £ f . = Z ^ = Z di 

3 . a 1 ^ 

U >0, j = 1 , . . . ,f.vi-K*-l 


(4.16) 


These constraints do not however pose any restriction 
on the polarity of the coefficients of the deconvolution 
filter, i.e. they will contain both positive and negative 
values. The only restriction on deconvolution filter is 
that- its order K should be small such that 


K « M 

where M is the length of the observed data. It has been 
observed that a deconvolution filter of order 3 to 5 work 
well in this context, for the data of different lengths. The 
constrained optimization of the objective function was 
achieved by using a NAG library subroutine EO^UAF, 



CHAPTER 5 


. SmULATION AND RESULTS 

5,1 FRIEDEN’S MXLYiUM ENTROPY RE3TORA.TION ; 

^ ) One-Diniens ional Dat a; 

Frieden* s maximum entropy method, discussed in 
Chapter 3, has been simulated on DEC- 1090 computer system. 

The object image data of size N, consisting of 32 discrete 
quantized gray levels, is degraded by passing it through a 
diffraction blur filter (sine (x)), of constraint length L. 
The resultant degraded image is of size M = N-fL-1. Gaussian 
noise with zero mean and standard deviation o is generated 
and added to the degraded image. 

From equations (2,15) and (2,16) we get a system of 
M+1 non-linear equations, which have been solved by the 
generalized Kewton-Raphson iterative method. In this, we 
start with for m=l,-. *,,M, and computed from 

equation (2.15), NoW the values of m=l, 2, . . , ,M+1 , are 
changed in a direction which satisfy all the M+1 non-linear 
equations given by (2.15) and (2.16), From the Nev^/ton-Raphson 
method, we obtain the new set of values previous 

set of values as follows: 
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r 

^1 



'^■2 

• 


' •' 1 

• 

• 

'^'M+1 

J 

new - 

. : 1 

\+iJ 


old 


[DEBAF]"*^ 


F 


where F is a vector of elements F^^ i=l,,,.,M+l 
by, from eqn. (2,16) and eqn, (2,15) respectively 


N 


A 


n=i 


for m = 1,2,, . . ,.M 


and 


^M+1 


Z O(x^) r Iq 
n=l 


The matrix [DEE5AF] is a (M+l) x (U+l) Jacobian of 


respect to the sst {K} f and is given by 

^F, 


[ DEBAF] 


dFi 6F^ 

5X- 575 


d F, , , , 6 F. , , 

M+l M+l 

07” * * • 


1 




^^M+l 


M*}“ 1 


(M+l) X 


(5.1) 


and given 


(5.2) 


(5.3) 

F with 


(5.4) 


(M+l ) 
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The iterative scheme of eqn. (5,1) can also be 
expressed as 

a = - [DEBAF]-^ F 

Where, 

^1 

The above equation can also be modelled as 

[DEBAF] a = "F (5.5) 

This expression represents a system of (M+l) linear equations 
in (M+l) variables , and can be easily solved with the help of 
NAG Library routines. In this case, FO^ARF has been used for 
solving (5,5). From the solution of (5,5) we can find the 
new set of values for {■*.] by 

(5.6) 

old 

A satisfactory set of values of , m=l ,2, . . ,M+1, is 

obtained after 10 to 40 iterations,. However, in this 
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implementation, about 15 iterations have been found to be 
suff icient. 

Selection of constants S and B; 

As discussed earlier, (2.8) and (2.9), the constant B 
should be such that it makes all the noise sample values 
positive. Since, it is not possible to have a-priori infor- 
mation of each and every noise sample, nor it is possible to 
find them deterministically, the most convenient valCie of B 
is taken as -2xo, where cr is the deviation of the white 
Gaussian noise process. It has been found that the value of 
B contributes largely towards the convergence of the algorithm. 
If the value of B is chosen larger than recjuired, the maximum 
value of residual, i.e, F^, after same number of iterations, 
also becomes large. On the other hand, if B is too small than 
desired, the process does not converge at all due to negative 
values of . Thus, the optimum value of B has to be found 

by trial and error. In the present simulation, however, since 
noise variance is known, the value of 3 has been kept equal 
to 2 X a and has been found to be quite suitable. 

The relative smoothing of noise and the recovered 
object image data is controlled by the constant S (2.10). 

For large values of it has been observed that the noiso is 
smoothed out more and hence in the restoration the edges of the 
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object become sharper. For smaller values of j the noise 
is more abrupt and the edges in the object are smoothened 
out. Generally, a value of - 20 lias been found suitable. 

Results ; 

The simulation results fox fevj impulsive objects as 
well as for randomly stepped objects have been obtained on 
DEC- 1090 computer system. The restorations for various 
objects for different values of o in N(0,o), and ^ has been 
obtained and attached as Appendix . B. Tlie results so obtained 

tally with the results published by B.R. Frieden in his 
paper [1 ]. For restorations obtained by restricting the 
number of iterations to 15, it takes about 9 second of computer 
time on DEC-1090 system for N=16 and L=17. For N=24 and L=25, 
the CPU time requirement is about 20 seconds. 

It can be seen from the plots of degraded images and 
their restorations, that, if the degradation of the object 
data sequence is considerable, as is the case here, then the 
smaller peaks get smoothered out, and also if the interval 
between the two successive pulses is small, then they get 
superimposed in the restoration. Also, as the size of input 
data decreases, the restoration becomes increasingly smooth. 
For this, two objects of same shape but different sizes, i.e. 
N=16 and 24, were taken. The restoration for Nrslb was much 
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more smooth as compared to N=24. Furthermore, the restora- 
tions obtained for impulsive objects have shown better 
resolution as compared to the randomly stepped objects, 

b) Two- Dimens ional Data ; 

The tv;o diinens ional object image is represented by a 
matrix of order N, which is degraded by a two dimensional 
diffraction blur filter, of constraint size LxL, with 

o o 

seperable impulse response (sine (x) sine (y)). The resultant 

degraded image of size IvlxM is further injected v^/ith Gaussian 

noise of zero mean and deviation a . For the purpose of 

processing, the two dimensional data is converted into 

one-dimensional data by vector-space representation. Thus, 

from (2.15) and (2.16), where N is replaced by N and M is 
2 2 

replaced by M get a system of M +1 non-linear equations, 
which are solved as before. 

Results ; 

For simulation purposes, an object image of size 8x8 
was taken, v^hich represented two bright spots against an 
uniform background. This was degraded by a filter of const- 
raint size 5x5; and noise N(0,o) was added. The resultant 
array of size 12x12 was processed to restore the original 
object image. The value of jf^ was chosen ’bo be 20. The 
convergence vjas obtained within 10 iteration to the order of 
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-4 

10 • The CPU time required v^as 3 minutes. The arrays 

involved for processing this I2xl2 blurred image are 

of size 145x145, This, as Vi/ell as the iterative nature of 
the scheme accounts for the excessively large CPU time 
requirement. For the sake of comparison, the computa- 
tional time requirement for the processing of a degraded 
image of size I6xl6 is of the order of 40 minutes. 

5.2 ONE*. DIMENSIONAL SEARCH ALGOftITH^A FOR MAXIMUM E^^TRDPY 
RESTORATION : 

One- Dimensional Data : 

The computer implementation of this method involves 
pre-selecting a value of p, and then finding a suitable value 
of A which maximizes the entropy measure (3.8), To achieve 
meaningful results, p, must satisfy constraint objective 
function (3,7). Based on the quadratic form (3.20), the values 
a,b and c are computed for computing the set V, (3,20a) and 
(3.2Db), and the set VJ, (3,^c) and (3,20d) , Now, the value 
of p, is chosen from VflW, This initiates the algorithm for 
the restoration of the degraded image in the presence of 
noise, as explained in flow chart of Fig. 3,1, 

The object imago data of length N=l6 is passed 
through a diff raction blur filter of constraint length L. 

A Gaussian noise with zero mean and deviation a is added. 
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This degraded image data sequence of length M=32 is 
processed to restore the original object data . 

Results ; 

The processing for various types of.-objects indicate 
that the convergence rate for impulsive objects is very slow 
and also the restorations so obtained are smoother than 
desired. However, for randomly stepped objects as well as 
for objects with gradual variation in intensities, the 
restoration results are fairly good. Also, the convergence 
is much quicker for such objects. The computational time 
requirement is considerably less for tliis algorithm in 
comparison to Frieden's aigoritlim. The increase in the size 
of the data being processed does not increase the computational 
time requirement in the same proportions as in the case of 
Frieden’s scheme. However, the resolution achieved is 
inferior to that of Friedon' s algorithm. The results obtained 
are attached as Appendix B, 

Tivo- Dimensional Data ; 

As before, the two dimensional object data is degraded 
by passing it through a two-dimensional differaction blur 
filter and subsequently Gaussian noise is added. For the 
purpose of processing, the tv^/o-dimensional data is converted 
to one-dimensional form. 
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Re sul ts : 

An arbitrary object image data of size 8x8 was degraded 
by a 9x9 diffraction blur filter . The resultant I6xl6 image 
array was processed to restore the object image data. The 
computational time requirement was of the order of 1 minute 
for an object having smooth variation of intensities. However, 
for impulsive objects the computational time requirement is 
much more due to the slow convergence, 

b.3 MINIMUM HdTROPY DECONVOLUTION: 


One-dimens ional object data sequences of size N, which 
are impulsive and hence are intutively known 'to have non- 
Gaussian probability distribution, are degraded by diffraction 
blur filter as before. Our aim is to design a deconvolution 
filter which, when convolved with the degraded image data 
restores the object data, minimizes _ the entropy or laaxi- 
mizeo the order in the estimated object data (4*13), The 
objective function used for the restoration is varirnax norm 
defined by (4,14) i.e. 


Ov = 


= S 


j • 3 


where fj 


is the estimated object data and is given by 


A M+K-l 
f . = S d. . 
^ i=l 


< 3 , 
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in which is the degraded image data for 

and gj^ is the impulse response coefficients of the decon- 

volution filter of order K* Substituting f . in the expression 

W 

for Oy, we get 

M+K~*l ^ M+K“l o o 

0„ = 2 ( 2 d , gJ7( 2 ( 2 d . g y 

3 3 i^l ^ 

This objective function is maximized under the const- 
raints 

M+K-1 ^ M N 

I f. = S d. = s 

j=l J i=l ^ n=l 

A 

and f ^ 0 » 

for designing the deconvolution filter coefficients g^^, 
k = 1,..,K. The restored data is then computed from (4,13) « 

A 

The constraint f . \> 0 ensures that the pixel values of the 
restored object data are positive. The other constraint 
preserves the amplitudes. 

Results ; 

The simulation results show that for very badly 
degraded images , the restoration of the original object data 
is excessively smooth. However, for the images which are 
degraded to lesser extents, by choosing the degrading filter 
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of smaller constraint lengths, it is possible to recover the 
original data by minimum entropy deconvolution. It has been 
observed from the restoration results that the deconvolution 
filters of smaller order, 3 to 5, give much better resolution 
as compared to the deconvolution filters of higher orders. 

This is understandable because the deconvolution filter tends 
to increase the spread of the estimated data, and thus the 
deconvolution filter of smaller order produces better resolution. 

Since the solution to the deconvolution problem is 
obtained iteratively and also the knowledge of point spread 
function of of the imaging system is completely unknown, the 
convergence of this scheme is siovv* To process a one-dimen- 
sional array of length 32, the CPU time requirement on the 
DEC 1090 system was of the order of 50 seconds. The various 
results so obtained are at'tached at Appendix B. 



CHAPTER 6 


CONCLUSIONS 

The ideas of maximum entropy method, for the resto~ 
ration of images, are intrigining since they suggest that 
one can in principle remove the presence of artifacts, reduce 
the apparent noise and give extra resolution. Maximum 
entropy restoration technique is more fundamental than the 
other techniques and deals with the problem in a very general 
way. What the maximum entropy method sets out to do is to 
seek the most likely answer to a particular question that 
can be obtained from a given body of scientific data; Entropy 
maximization is a kind of insurance policy that protects 
the restoration against spurious de'tails for which their is 
no evidence in the observed data. 

Amongst the tvjo maximum entropy approaches considered 
in Chapter 2 and 3 it has been observed that the Frieden’ s 
maximum entropy method gives better results in terms of the 
resolution achieved as well as the rate of convergence of the 
algorithm, whereas the one-dimensional search algorithm is 
computationally faster for larger arrays. The efficiency 
of the Frieden’ s algorithm lies in the fact that it seeks 
the minimization of the error in each pixel of the restored 
image, whereas the one-dimensional search approach minimizes. 
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the average error in the restored image. Also, the conver- 
gence rate for the impulsive data images is exceptionally 
slow in case of one-dimens ional search algorithm based on 
the solution of initial value differential equations. The 
simulation results indicate that the one-dimensional search 
algorithm is not efficient enough "to offset the disadvantage 
of excessive computational time requirement in the case of 
(M+l) dimensional search algorithm. 

The ideas of minimum entropy deconvolution apparently 
dqes' not offer any promise in the field of image processing. 
But, its extensive application in the field of seismic data 
processing does indicate its applicability in the field of 
image restoration for a class of image signals which are 
impulsive and have essentially a non-Gaussian probability 
distribution. The principle of minimum entropy deconvolution, 
for Such a class of image signals, has been applied for the 
restoration of the object images from diffraction blurred 
images. In the absence of any a-prioxi information about 
the point spread function of the imaging system, the results 
obtained are quite encouraging. However, if the degradation 
suffered by the object image is very profound, as well as if 
the noise level is high, the minimum entropy deconvolution 
method fails to recover the object image. 
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APPBiDlX A 


DIGITAL IMAGE CHARACTERISATION 

In digital image processing system one usually deals 
with arrays of numbers obtained by spatially sampling points 
of a physical image. After processing another array of 
numbers is produced and these numbers are then used to 
reconstruct a continuous image for viewing. Image samples 
normally represent some physical measurements of a continuous 
image field. For example, measurements of the image intensity 
or photographic density. 


For the purpose of analysis, it is convenient to convert 
the image matrix to vector form by column or row scanning 
and then stringing the elements together in a long vector. 
Consider an image matrix F such that • 


F(l,l) F(2,l) F(M^,1) 

F(2,l) • 


(A,l) 


F(Ni,l) 


F(I\,N2) 


NiXN^ 


An equivalent scanning operation can be expressed in 
quantitative form by the use of an N^xl operational vector 
and Nj|^N 2 XNj^ matrix defined as 



A. 2 
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^2 ’ 
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(A.2) 


v/here each element of represents a matrix. 

Then the vector representation of the data matrix F is given 
by the stacking operation 


N 


f = 


2 
2 

n=l 




(A.3) 


In essence, the vector extracts the nth column 
from F and the matrix places this column into the nth 
segment of the vector f. Thus f contains the column scanned 
elements of F. The inverse relation of converting the vector 
f into matrix form is obtained from 


5=2^ f (A. 4) 

n=l ^ ^ 

VJith the matrix to vector operator of eq, (A, 3) and the 
vector to matrix operator of eq,(A,4), it is easy to 
convert betv^een vector and matrix representation. The 



advantages of dealing with images . in vector form are a more 
compact notation and formulation of a tedious two-dimensional 
problem into a one- dimensional problem* Thus, the results 
for one- dimensional signal processing applications can be 
easily applied to two-dimensional case, 

FIInIITE area superposition operator 

Considering the discrete superposition of a spatially 
truncated data array F(nj^,n2) for nj^,n2 = 1,2,.,.N v;ith a 
spa daily truncated impulse response operator H( if 2^ m^,mp 
for 2 ~ In general., the impulse response 

array changes form for each point (mj^,m2) in the processed 
array Q, The finite area superposition operation is defined 
as 

^1 "^2 

Q(m^,in2) =2 S F(nj^,n2)H(m^-nj^4-l,m2-n2+l; 

n^=l n2=l 

(A. 5) 

for mj^,m2 = 1;2,...M, where H and F are assumed to be zero 
outside their range of indices. From the indices of the 
impulse response array at its external position indicates 
that M = N+L-l, and hence the processed array Q is of larger 
dimension than the data array F, 

In vector-space notation, the finite superposition 


operation can be written as 



q = Df 


A. 4 

(A. 7) 

2 2 

where D is M xN matrix containing the elements of the impulse 
response. It is convenient, to partition the superposition 
operator matrix D into submatrices of dimension MxfJ, such 
that 



h,i ° 

0 



"^2,1 °2,2 

• • 

• * 

• 

• 


D = 

« • 

1 °L,1 ^,2 

i 

m 

' .(A. 8) 

i 



• 


1 

0 0 

V • 

* 



• • 

• • 

0 0 

L. 

w 

,N 



The general nonzero term of D is then given by 


V'/here 1 < n. < N and n. < m. ^ n,.+L“l. Thus, D is highly 
structured and quite sparse, with the centre band of 
submatrices containing stripes of zero elements. 


If the impulse response is position invariant, then 
the structure of D does not explicitly depend on the output 



rray coordinate Also, 


Ha2,n2 ~ Hn2+l,n2+l 


(A. 10) 


Thus, the columns of D are shifted versions of the first 
column. Under these conditions, D is known as finite area 
convolution operator. 


For a specially invariant and seperable impulse 


response 


H = hpj 


(A. 11) 


where hp^ and hQ are column vectors representing row and 
column impulse responses, then 


D = X Djj 


(A. 12) 


The matrices Do and are MxN matrices of the form 

n. C 


°R = 


hp^(l) 0 
hj^(2) hp^(i) 
hp(2) 


jhR(L) 


hpa) 


hR(l)| 


hR(2)| 


hR(L)J 


0 


MxN 



The two-dimensional convolution operation may then 
be computed by sequential row and column one-dimenSional 
convolutions. Thus 

Q D^F (A, 13) 

In vector form, finite area superposition or convolu- 
tion opera ter requires nV operations if the zero multi- 
plication of D are avoided. The seporable operator can be 
computed with only NL(M+N) operations. 



APPENDIX B 


The following pages B-2 to B-123 are divided into 
three subheads pertaining to three different image resto- 
ration techniques under perview. Each subhead constain 
the FDRTRAN-10 listings of the simulation programs followed 
by the plots of the simulation results. The various 
programs and the results contain sufficient comments for 
easy assimilation. 
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PROGRAM BARK. t OK ROk uEGPrtDIi'jG THf. -jBjECT iMAuF OhTk, 

TRIG PRuGRA'i JEGRAOhS TijF OBGLCi IMAGi- ukjh P.i pac.StMG 
TT THHO' A rUFfRATinrt dLUP Ki f. IFk (oT .-jC (. X j + J . jt'RG FTLI'GR 
COFFFS aRG 'jnRMAT.lGKn SiJCh THAT TnPil' SU^ TG gO'JSL J'O t. 


program BbiJR.t'OH 

I WThGG R t C 1 1 o ) , I ;'iG C B i A , LX , H a 

REAL LT;1T.,r,jL'lK,dR(171 ,OK(G7 , lo) , ALRu A , SUM , X O-G, ( ,1 ? J 

MX = i6 

h'A-ll 

MX=wX+GX-i 

0 P E W C H ;i T 1’ = 2 1 , n R V I C E = ' H A K ' , r T L E = ' p T C . 0 A T ' ) 

RE AtJ ( 2 i . * ) , C PIC C X ) , T = 1, , 'Gv 1 

LlMu^-B.lAlGG 

LiMRs J . i 'USy 

SuM=f<.<- 

nu iR T = 1,T,X 

ALPHA=LI 'i;jtl-<'CL.i;fU<-LTHL) / LX 
H K ( i ) = f o T < I ( A L P H A 1 / A t, P t i A J t + ? 

.SUH = Su:itTKCl') 

COHriiiUE 
DU 20 T=t,Lx 

COU'i'TiP'G 
DU **<' 1 = 1 , MX 
DU pR 1=3, ,Na 
DK( i f J )=U .0 
CONTIRHi^ 

CONVIuHC 
DU UL >1 = 1 f H a 
DQ 7R T=1,UX 
DH(i+>J-l , J) = Hii(i') 

CO'J‘i'lL'’E 
r U I'l 'l’ T I'H ' I*. 

WRITE ('45 ,Ry) , ( U)R ( i , J ) , >3 = 1 , Ha) , T = 1 , HO 
FuRiiAT(lLtE10,L» lA ) J 
DU i('U i = i»HX 
SUM=l>.(‘ 

Du 11 L .\=l#i<X 
SuM=^SiJ■'f OR( I ,K) *P LC (.K3 
CUfiTTijnc 

TMGtl J=XElXC.SU?0 
XiMUCl )=S>1M 
COMTTililb 

' WRTTEUiOO^R) (IMG(I),T = 1 ,MX) 

WRI TP_^ 5u , '’* ) ( A I !’iG (. I ) f 1 — I r i‘‘X J 

FuRi-iATCieCTy ,1X) ) 

srop 



TliK Pu'U/^a'5 adds 'fllUTti; ijAuS;i.U^ ^CCfSiniViA 

Tu THh iihu^HEO i'lyvCti ijATA, 


PunbRA'”^ 'iul tiK . H uR 

IfjTt:Gr:P T. iG t 32 ) , MX , ijX , U , 1 SUM , P iC (. 1 b ) 

RtAL .■>T(,!lA,iiRAG,S;jM,xr*U(32J ,S^iK 

riX=i7 

MX=nX+LX-i 

npEu (HuTT=2l ,ni':v,iC£='nsK * ,ktlf='p:[c.ijai ' j 
0 p E 1 1 t Mi 1 1 X’s 2 2 , n £'ri c £= ' D S K ' , £ T, uE = ' £ M j< 5 U , A i ' j 
REAUC21 ,*J CP.t,C( t j , 1 = 1 ,ijXj 
Rt;AM(2?, CXlMdCn ,T = 1 ,’1A) 

ToMii=u 

!>u IM T = 1,''x 

ruN'i'CtlUh 
SiGilA=M, 1 u 

ISiIh=XSu i/ 'lX 

s R = I ii ♦ A r. U G 1 r; (. I j; M m / S T U A * * 2 ) 
nu 2 '* Tsi,'';! 

SU M-X 1 ‘U., ( I ) tGOS UOF ( i-iF , S iGi iA J 

ItiGCI J = if’’iXlSij’!j 
XX2(j(i3 = .S.J ■( 

CU’lTIiJDt;, 

WRlTf:(''>l , JO) (ThGCI) , 1 = 1 ,f’X) 

WRITE t S 1 , ^< ) ( Al bG t. T ) , .1 = 1 , hX ) 

WRITE (.57 , /7) 

FURfiATCiOA, 'SAP IS =*,EiO.'^) 
FuR'-*AT(lOap,lX) ) 

SIHR 

end 



0 m» 5 


Pk'HaRrt'l FLIH FRiKuF.iM'i, HAAT.^th*! t^TJl'ktjPi I -J^'^ 

TtlTS PRGGK'Vi-i SOLVtS A SYSlFh uF A'u;>!-I.iWt;AK hOuATIuA'o in 
RESrOHE TiiE nuIijM^L nuJLC'i IMAGE OMA FkO’, GEGRADEn 

TMA'jE H/iTA, 


PRHGRmM MK.FOR 

REAE B,P(i’5) ,uFC 33,3'^ J ,Lf33) , SU '2 , GI’a 3 ,Rij , ATG'’'A , -’.AX 

, G ( i2 , 1 o ) , SUM U 6 i , Tf; , [MG ( 3 7 ) 

Integer .‘.ug c ( 32 , 2 j , .irujF i , ,-a , ?xx , ur , -'m , 12 , -'3 , u , 1 c ( 1 6 j 

C (,l M'3n M S f o r li T f ^ # *.5 i* # E f *3 n r* # .■ ' 1 A X / i. G r ('-) # 'XM f f /V f ij X » M ^ X 

SlGwA=:M.2 

B = 2.U4,SXG-iA 

MX"1G 

1(X = 17 

+riA“'i 

MOnETslG 

R0=-2A.0 

DPEn ( .ji (i T=:2 1 , i;E V TC E= ' liSK ' , F n,h= ' V-'0R5 \ .MAT * ■) 

PEAG C 2 1 , * ) ( TiiG ( T j , 1= I , mX ) 

npE;'i ( UA! t T = 2 2 , uF M ICr= ' OSK ' , FiT.E= ' FjRaS .H^T ' ) 

neAD(22,tj( CSCI,vl^ ,J=i ri'fX) ,i = l,iiXj 

S Tb THl CO^a^LOTifA, aAXRxX CF THE IMAGinG SYSTEM, 

DU l'> I = t,l-.X 
siHTCia ) = r 
siWi.'Ci;,2)=t 
ruPTiuML 

SXN'i’iuX,! J=LX-1 
DU ^0 I-liX+i,fiX 
RIfri'Ci,l)=MX-l + l 
SI Ml' ( X ^ 2 ) = T "•IjX+ 1 
CUfiTIUMl:: 

TU=U.G 

DU 30 .T=t,1X 

TO = l'> + T'lGi.J)*l.C- 

rofJTIii'lE 

DO 4(* T = 1 , /ix 

T.Ct J=O.G 

ru'Fi'TiJUt. 

0 1 M A+ i ) =- 1 - AOuG C T 0 ♦ 1 , U / iJX ) 

DU b'' 1=1 , "iX 

SUf'UT J=0Xi'>(-l-?aU+ 1 )) 

ru'iTiuno 

CAl.O GEIF 

Ckhk^ GEXDES 

no UD 1 = 1 . ,'ionFT 

CALIj gexla ^ 

CAlUi G EXSUM 
CADu GETF 



ou IH> i-lf-.i'i 
PlCa)=iH-'iXt.SUMCTJ ) 
cuirriH'.’L 

WHT'i:K(4H a'SO) CP ICCIJ , 1 = 1 ,1' xj 

FOHiiRTf K>tlb,ixj ) 

srnt> 

FtjD 


surmnuTi'iii g»;:tf iHt HF-^TuMnL. ne’ THh cn^^SiP/iiiv Ti, 

AFTC.R JLTLRA'rir);i.K(i1,T=l..,.|Xft,.Sri1.jT,u HP N.fcAKr,x a'-hH 

AFTLR TiiP SpPCTFTGn 'ju. Ok ITKPAT'iniX, 


suru<nu'Fi‘it-. i.Ki'F 

P L A U S f ) -i 1 , S , T'" , D H , u , 6 : 3 ' I , M A a , i . , K F' , T i; , i. ^ 1 
IiJT'k.GFR SIH LMjX, ,K,Ka 

COM.'IOim ;> ( i? , t o ) , S I II J- f .j2 , ? > , K f 3 4 J , ,jF ( 3 3 , i 3 j , j, f 3 3 J , b'i I f I H j 

,HAX,Ifiu( 32 ) ,!X, uX , Fa, iO 

MaX = 0 Ji 

DU iR T=1,U 

Ftl J=P.'- 

DD 2 (' f 1 =l ,Siri'(i,n 
P c T ) =F C J.) +.S a , 3) ’*=^11/. ( J ) 

J=:J + 1 
COM'i’lP'lF) 

F c T ) =F C i ) +f:XP ( - 1 “ii ( X ) /Pij ) -P« ThF, C X ) 

I j ) Guto lo 

MAXs:F(Xj 

COn'i’TiiilL 

FtMA+i )=o.e 

DU 3 '' 1 = 1 , 3 X 

FtMA + l )=Ft 4 a + 1 ) +SU'UT) 

CUWJ.'Ti4Mh 

FC'U + l )='!' ( UFl ) - I'i 
T F ( H A A , uP . F" t A+ i ) J UnTO '1 f 
MAX=Ft'U+ i ■) 

CUiriTi.ni'.; 

Rt^TlJlFi 

FUD 


SuPnnuTi ji’., GFXnc.P CU’iPUXFb XHk, ..lACHttlA'^ in- vEl'^jP F HiTn 
PEF> HF-C X i'u 4 HE i*Ai‘iPA'o. 


RUPHO^i’fiME UF'IDEH 
X u X E G 1.4 H G X !‘i X , 1 ) X , I'i A f i 1 X 



'3-7 


, !»1 A X. f I M ij C 3 2 } ^ b , KO f I j X » ij X , ^’X » i 0 

DO iO 1=1, MX 

DO 2C J=1,MX 

DFCi,J)=0,C 

K=SIMTCj ,2) 

DU if- '1=1 ,Sl!lX(0,13 

DF C 1 , J ) =DK ( i , 0 ) -S C .1 , k 3 *5 ( 1 , K 3 ^SU’’ (.K ) 
KsKtl 

comtiwue 

IF (I. of:. 0 3 onxo 20 

DF' ( I , J ) -DF ( I , J ) - C FXP r -1 -T c T J /HH ) ) /Pj 
cof’ttwDj: 

CUDI-XnOK 

no 4 D I=I,!AA 
Dt (i,ilX+l)= 0 .'J 
J=S1MTCI ,2) 

no t>(> K=i ,si 'll ( 1 , n 

DF ( i , OX 1 1 3 = 0 F 1 1 , ? 1 4 + i ) -r, c T , J ) '♦ 5 D. 1 ( a ) 

,T = ,1 + 1 
CUfJXlH'll:. 
ruH'I'TiMfJr- 
DO OD T. = 1,'AX 

DF(rtX + l ,T J = (v.O 

J=SIMXCI,2) 

DU 7f' K = 1 ,S.L'-i'(l,i3 

DF (wX + 1 , T J = UFC?'''X + i , I ) -FhT , JJ *SD -if J) 

J=iT + l 

COMi'iHUl. 

CUflTTNUU 

DF (mX + 1 ,MX+1 1=-Ff nX + 1 ) +j.O 

RclXuRh 

F.'jD 


ThF SjDH'JuTi'Jt-. vjFTSilM CUN''FD'tF.b c-STloAiFu OPjECT DATA 


AfTi::P FACW X'^FHATino. 


SUnKDUTiMtX OKI'S UM 
I T li G L R s ,i r-i , 1 : , 0 X , !'U , r i X 

R ijj A O S 0 1 , S , F , D F , ij , S D '"i / • ’* A X , i3 , R D # T u , .1, ''1 0 

rOMiin.i S(i2, lb) ,5i;pi’( 32,2) ,F(33) ,DK(.33,33J ,uC33) ,S'J»i( i6| , 

M A X , I OG 1 1 2 ) , D , n 0 , ti X , 0 A , i-' X , 1 0 

nu ID r=i, , '>x 

SIJ'HI)=:D.0 

J=T 

no iiD K=l,fjX 

SiJII F T j =StT i ( X ) -L ( J ) 1=0 r J , i 1 

J=J + 1 

ClJM'i’TiJDL 

So '-U T ) =0U 1 (I ) - 1 -0 C X+ 1. ) 

SU'HI)=nXRCSi!Hf I ) ) 

CunTIuDO 



C,5t’'i;LA'’ StiljVF^ij a ^>VbT>'Jf* Of f-lX+l I* A e. P i i 
cO'-iFiirr: Tsip: .'Iajx' sct iw 


SuHKFUjTi'ff; viK'rL/vM 
I w T t G t; S t n J , u X / ■■! A , I -i X 

R A* A ( j iSGr! T f S f P f PR* / ]j ^ n ■'] » M A X f t5 / HR f IMj f A 3 j # i. A' ^ S iJ 2 ^ X u 
CUMMfVj r> c 3 , I b 3 , S If'X ( 3^ , 2 j , P ( 3 3 ) , t)P' t 3 3 , 3 3 ) , b C 3 3 j , uH, i ( 1 6 j 
r HAA r ,l'R, ( i2 ) , b , I'b'f , ‘-a » bX , '*X , I 0 

neAL XRf (33,33), XP(33j,ivKSPCPt3i3 

rwTi^Gt'R IA,N,iPAXb 

DO 20 U I = i,;"iXtl 

nu 21u vl = i,'iX + l 

XuFtI,J)=')f’CT,JJ 

CuM'i'TuOh 

CONTIb'U: 

DO 220 i = l,KiX+l 
XF(i)=:-l*PCX) 

CoDt’lonu 
TA=-)X + .i 
IFAiLsO 
M=MX+1 

Chh'u F f “i A R, !>' ( X D F , I A , X P , •'I , n i, , W K .3 P C P , T P A 1 0 ) 
IF(iFATb.wP.b)STnt> 

DU 23o i = i,.iXH-1 
r.CT j=bCi )+Di.,( i) 

CUOTTOMO 

RtiTURi) 

F.WD 



PROGAm PI.LIT.FijP TU pT.UT ThE TwAGE O^^TaS. 


PROGRAM PuflT.EfJK 
TNTGGGR , ;iX , r.X, FTC ( 32 ) 

REAL WR(6i) 

common la , T.X , iX , ^ IC , ••■JH 

A?X = i6 

T.X=17 

MX=i'iX+LX-l 

npEU(UMlT=:21 ,L)EVTCE='DSX' ,FiT E='r^ir,nAT' ) 

PEAD(21 ,♦) (PICC 1) , 1 = 1 ,A'X) 

PRTI'^T to 

EO'’ilA'rClH ,/////// ,40 X, 'The ObJCCi j/'MCr/j 
CAT..L PPLOT 

OPEN (UM1T = 22 ,i)EVT:CE= 'DSK ' , FIT.E= 'FUR4^' .OA'^ ' 1 
READ(2?,4J CPICf i) ,T = 1 ,MX) 

PRINT 20 

■ FORoATf iOi, ///////, iax.'niFKRACi’Ti,:'''' 'HiOKPED Tr-iAu’^') 
CAIjL I plot 

OPEN ( UM1T = 2 3 , LEV TCE= ' l)SK ' , FIT E= ' POP 4 1 .OA'^ ' t 
PEADf23,*J('PIC(i') ,T = 1 ,fU) 

PRINT 3U 

FORMATCl 11 ,/////// ,2'')X/ 'DIFFRACiTO''’ Si,UkPCO T.iA.;,P 'On I’h 
iHITE GAOGSlA..f ijniSC 
CAT.L IPl.OT 

OPEN C UN 1 T=24 , DE V IC E= ' OSK ' , F l.T.i':= ' FuP 4R . DAT ' ) 
PEADr24,*)(PTC(I), 1=1 ,MX) 

PRINT AO 

FORMATdHl, ///////, 20X, 'PK.CTnHP;,) j ‘'AG l f‘"'l,TKio-,NG 
MAXIMUM FijTRnpY PRT jClPoF ' 1 
c.^hLi FPicrr 
STOP 
END 


SUPROUTIME PPI.OT 

integer NX, ox, MX, pic 
real wn 

r 0 M M 0 0 N X , L X , i i X , P 1 C t 3 2. ) , F j > f b 4 J 
CAI.L SCAOECi'iX) 
print E 

FURMATClIi ,////) 

no IP 1=1,32 

DU 20 J=t,f>4 
Wti(J)=' ' 

COM TIN nr: 
no 40 K=1,NX 
L=32-l 

TFCPTCIn) .NE.ulGQI!! 40 
WB(1*10='4' 

COMTINUE 

PRINT 3U , (GbC J) , J~1 ,G4) 
FOPNATCIH ,iOX, * ! ' ,5X,64Ai) 
CONTI fj’U-.: 




1 »,* 


FORMAT ( IH ,1 nx , 1 R X , ' 

— ) 

RETUR./J 

FWD 


R(iRaouTi?‘if-: iPi.nr 

integer I X , NX , , PIC 

real r/R 

COM H n u i j X , h X , ! . X , D i; r ( 3 2 ) , ! ‘ ti C u 4 j 
CALL SCALE Ci'!X) 

PhTNT 5 

FORHA'i’Ciil ,////) 
no 1C 1=1,12 
DU 20 7=t,64 
WBCJ)=' ' 

COMTTfHlE 
no -in K = t,«x 
T. = 32?I 

TFfRTC(X) .’SE.DGnm 40 
Wli(2*K) = ''*' 

ClJWTTNITE 

PRINT 30, C -laCJl , J=1 ,6t) 

FQRi'iAfflH ,iOX,'! ',SX,r->4At) 

cqnttnid-; 

PRTN'^ So 

FUH ‘i A T C 1 HO , 1 ox , ' I ' , / , 1 D X , ' 


RKTUfPN 

FND 
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integer M 0 , N X , L X , 1 X . P I C , »■’ I N , X 

real NP,S1JM 

CUMODN X , LX , OX , P tC C 3 2 4 , Wb f p4 ) 
MAX=0 
M j;N=0 

DO 10 T = t/IU 

Tb’CPICi n .GE.hI.OGOTO 20 
f4lt'5=PIC(T J 

TFCPTC ( I ) .LE.ilAXlGnrO 30 
maX=R1CCI) 

CONTI NOE 
COMTTW"E 

TFCMTW.GF.OGUTU 'in 
no 50 1=1, MU 
DrC(.T) = PIC(il-MIO 
COOTTNUt-: 

max=max-m:in 

TF(MAX.LI-'.31)l,nTn hO 

nu 7 C 1 = 1 , MO 
y;tjM=PiC (T ) *31 .0/MAX 
RIC CT)=IRIXCSU'4) 

COMTIUHE 

CUNTXiO'E 

RETURN 

FND 
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PRnuRA'’ -IbftaO.FUR 

l>lMh.»f6(UM !)l(( 12,1) A), b( 144.04 45) 

niMLMbl.jM Dc.HAKll45,l45),|,n.MnftCl4!>),RinAfb4),XT*'iGH2,l2J 
f^clALi ti.lf ,^)FriAK, 1 . A ■■iDA , S lGi^iA , aX , , I)K , tjllr'i , o , iMu 
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C IrillTAM^AM'lO.i U**' FSTlM/iTEP PaTA 

nu 11 Tal,‘lX 

.Su‘UT)»i.XP(-l-LAMUMl'>A + i) ) 

11 CU-lTT;4:Ft: 

r FURrtfJuAlIOfJ OF rMii, CO»VuT,UT10h- NM'OIX S 

nu 2f lal,’4A 

KfaX 

CaI^u rta*(KT,Mi,H2,M3) 

nu ic T8l,64 
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SOM i T J sS=» I ( 1 ) - 1 -L A -ID A ( 1 45 ; 

SUM(t)ai!,Xi;>(&:in(l) ) 

to coMVTjni!; 


C .SO*»Rno'l'i.Mt. ufirliMM COMPIIIFS THIS liF* iIAIjOK of IiA^^UA'S Ml 

C SOTiVImJ; a S^T or' t41i iiTliFAP Fynii’^in.iS, 

c 

SiiHKnuTlNc: OKILAM 

RisAlj l)rj(t «5},UKbPCF.(U5) ,S,K,uFuAF ,iiAf>inAf filli<i,:lAX»10,Ti4r 
n&Aii Xr(l'*5),M,Kn 

r v,H,l®‘Al U,8T«T 

CuMrtn.l ^SiiOh (o4 J ,h ( i4'j) ,ri’jOAF( 1 45,145), iiAftOA 

1 (146) ,80 '41^,1140.(144) ,i),Hri,r.(i4t,r>4)»sx?4;ir)44, j),Tt 

nil 10 Tal,l45 
XF(I)a-K(l) 
to COR'i’Tri'lii 

lAai4j 

TFAil.an 

f!al4l5 

CATib KO*A«F(ni£l)AF,TA»XF,»',niir«»FsrCF,TF.M(.) 

TF(iFAtu.*4E.O) 6TijP 

no 30 1x1,145 

r.AHOA (T)ai,AHnA(l)4-nij(13 

20 

RISlTURi't 

FM!) 



C? 




H'liiii KI.M!*h rHfi tLfcfc'ISMT.'S OF HHS KHn^FAbSc;') Ab 
T^PoT Prtt^/iiltiTiB^PF tUft CnhVul JT11“ MArKTA S,IT AijRu Fi«ns 
TmK NOI ZdBij p,l,i;rth»JXS Ttt KACH TU aVOTk h'liiTlPijTCAXlO-Jtf 
-UTti /iRlCl KLiRwEuTS. 


SURKPUTlBii Bflft(Rr?jir!u,ftl ,**‘i ,hM 
T Ri f Ml f i'll? » J 

f'cl Aij «JR f fiKOirf tit OR(J Af fiSM^i#D#ilAX#XPf TiiRj T iAP I) A / Pu 
m liRU2,a5,SrtOrt(o4),t (14t>1 ,ni;;RAF(l4S,l45j,liArtnA 

i a4‘'0,a'>:i(»>4 ),l•IAA,lMl;(l44^,R,Ro,<>(l*'*#»>4)..ST.<T^l^J,3J,.lft 

Pu '*Pw xaifu4 
ShnM(i)si).0 
400 rvj’iiTMrii., 

fflatHj M.jn-u/12 
K;4SKn<<=^ii«l?*Ki 
K1S\1 +1 

TK(ivl .Rv,!>) GUTIJ 410 
TF(Al .RO.O.ijR.^.nu.b) GUTU 42(i 
TF(^1 .r.i.bl GOTU 480 
M13b«(Ki~b) 

XaKl-b 

(Jgflj 440 

410 Mlsb 

XaO 

CtjTiJ 4l0 
420 Miab 

XaKl-b 
GuTU 440 
430 

VaO 

440 TF(K2,l5y.5)GntO 430 

TF(K2,Ry.b.UR,7.ni<,b)OPTO 459 
Tf (K?.l»1.3)b'3"iO 400 
M^a3-iK2"c») 

M38H*X*K2-9tl 
YxF2-b 
noTt^ 47 u 
455 «2ab 

wia»j4A*A2-Stl 
YaK2-b 
rtUTiJ 470 
450 M2a^ 

wjao*X4l 

YaO 

GOTO 47o 
450 M2aK2 

Mjadbx+i 

yaO 
.laMi 


470 







Msx+1 

T.2Su1*mi-i 

r.jai+i 

nii bCu ^sLt,L2 
Du dio u^u3#Ij4 
.<;KnM(J)snK(Kl,K)*l)K(K2fii) 
TsJfl 

510 rOWVTijIlK 

.Ta.TtO-'li 
500 CU-iVTuUii 

l>l::TUHij 
KiiO 
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THIS PQHSkAi PTjOTK.Fnit PLOTS TriR TuiO Dl Mfc,)iSTa^ AL PiCT'JttES, 
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100 


110 

120 

130 

140 
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160 


170 


Tf4TKr;ER*2 TA 

ni^ywSTlJ'! TA(ti»,64) ,T|)(«096) 

npFn('j»iiTa2l,l>FyTCEa'USK*,RlTil!;e*FOR32,r)AT*) 

NlSAt;(21,*),((lA(T,J),Jslf64jasl,b4J 

NXap4; HKS64; I<AW«1 ; TLsQ ; IHs.^ i r i^FGal ; I iRb42 

rAT.L USP(lA,7X,nYftiAH,Tli»lH,NiC6«T.G) 

rKTURJi 

Fan 


SURAH JT £\'fa: U8P ( A A, HXf I4Y , I.AW « Tl, , 1!? , MtlG ^LG ) 
riJ>4i‘tOif lAlA4,b4) 

Tr<TliGEP*2 TA 


TNTUGtriP*? THrb1,64) ,rj|bVC32) tOLAi>)K(5) 
TiUGirAI.*] M4R(12»,5),nHAyr3?rRj 
niM|jf48TirT TAr61f64) 

PATA GRAy/6**rt*,5**H','X','ii*,'X'.*tl' 

1 'F 'O', 'S', 'b*,'! 


+ ' 


2'f', 's', ' 
3'*', 4*'-', 5*' 'r'B 
43*' ','-',3*' ',4*'*', 'O', 
S24*' ','a','0','U',2!**' 

0NBi2,Q 

FLBlti 

FHBlM 
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9 
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IF((F»1-FL).GT,0.0) CO Tu lOO 


TbFL 

FLBfH 

PriBT 

RAN6EB(FH-FL+1)/Gri 

AA»(S0RT(KH)-SQRTIFL))/GP 

FE*(PH-Fl»)/AT,iJOlOry+l,0) 

TaAtiAXl (Rij,laO) 

SSB- f I . 0 /«■ J ) * AI.OG C F II /T ) 

00 16(i isl,32 

CD TO (11<;,12»j,i3c,l4(;),l>Ai^ 
FLFV*KIitCl**U*RA*^GFfO,5 
Cfj I'D tbO 

FLEy*(SJ>’T(Fti) + (T-1 J*AA)**2tO,5 
CO TO ise 

FLEVaFT,tFli*AUOCCFLnAT( T 0+0.5 
GO TO 150 

FLRy*FH*EAPC5o*(G i-I)J+0.5 

LBV(T)sFL«V 

C0OTTMIU5 

IPCnX.GY.61) i«Xs 64 
TFCi«y.GT.fa4) .<Ys64 
DO IRO XBl,WX 
DO 18u jBl,f<iY 
FLTbI 

DU I7u Kal,J2 

IF CIA(l,J).Gl!:.LFVCR)) KMaK 
CO!«TINDK 



IttCl f J)sKuT 
tbO CDMITMlIi*: 

VRTTECLG^l) 

IXanX 

TY»2»riY 

OU £«lflX 

no 1*>0 jst,5 
rlankcj jso 

1-90 COMTTiini^ 

no 200 KB2flY,2 

.7«K/2 

NtiamClfJ) 

TFCH^G.IilCt.O) MCa33-f<(! 
no 200 ii^lrS 
OI^'e:(K-l«b)a<>HAYC>v6,OJ 
r,r*e(K,b)aGKAYCfir;,b) 

TFCNG.NK.i?) uTiA»KCii)Bl 
200 CUMTTi4U|j 

tlKTTBiTiC.»2) 
no 210 L 

TF(HriA»fKCu).SO.O) Go ‘1*0 210 

white CtiG»3) lTjI*l|£(M,b),Wal,TY) 
210 COMTIMOfe: 

1 PORHAT(lHl) 

2 |P0R«AT(1‘< ) 

3 PORWAXd'l*'# lX»12?<nl) 

RB'OORN 

E«n 
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GflTROPY RESTORATIbB RX Sni,VtRG A SYSTEM OF lMl'riAi4 wAti!»E 

Program mer.fur sbarcubs nh^j^y. a «.F:i,L 
IS A« oNfi DlMEiJSlHAli SEArCh pRo«be.i.rHE 
M«I,!iImIuH BETjUNGS Tw AN OPEiJ INFTnll'iS T^TERtfAb , 
TS tOllND BY trial PQR QUICKER COMVEHCeMCE.FUR TdlS «EaSu!JI,T.|K 
VALGE Of Al] AmD the NQ.bF I'fRRATlUNS IS Ftl) TriPU'IGH A.l hCLTpT 
SIAYErBMT, THfc input UAIA TU THiS PROGRAM IS tHt C0.«VuIiuTjl()t4 

***I?i*-^^ IMAGING SYSTFr ANu THE r4niSE INFECTED DECiiAoEu 
niFFRACTION BLURRED lAAGF. DATA. 


PR0LRA1 RER.PJR 

DlMEMStUN A (64 f 641 fU(b4j«F(b4) , PF(o4 J , ATt64,o4 1 , A PaCo t »(i-4 ) 
DIMEHBlul X(64)fXX(o4J#AMAT(S0,.E) 

PEAL iiAMn,LA?1r*,rtO, TNCR, EPS, gP, SIGMA, SGrt,AA,.BrC»T 

TwTEGER N,m,MOIX,PIC(o4),L 

CbNhON/ARRAXS/A,D,F,Ff,ATA,AX 

CGMMOa / V AR/gP , [i AMU , L ABN 

CUMMnN/CQMST/t4,rf,L, SIGMA 

Nslo 

LBl? 

M«N+L-1 

S1GHAsD,25 

EPS»n,0i 

0PENCuNXTa23,UEVICE3*|)SK',FlLEa'F0R34,DAT*) 

npEM(uWITa25,uEVICE«*OSK*,riLEa*FOR36,DAT*) 

REAU(23,«) ((A(I,J),j8t,N},isl,.4) 

ACN,N) IS THE CONVOLUYlUN MATRIX oP TriE IMAGT:«G SXSfE4. 
R&AU(25,«) (0(I},181,M) 

D(M) IS THE DEGkAuEU IMaGE rTTH NOISE 4(0,0,SXGriAj . 

TbO.O 

DO 5 181, i . 

T«T+0(I) 

CUN'I'InHE 

AAaO.O 

RsO.C 

C*0,D 

no 10 j8i,M 

SUM80.0 

no 20 tsi,n 

SUMaSUMtAU,!) 

ruNXiNtlE 

BaRTSOMADCJ) 

SUM8SUM**2 

AA^AAfSUM 

C«C+0(JJ442 

continue 

AAsO;S4AA/SlGiiA 

Bs-S/STGMA 

C80.5«C/S1GnA 

DISC8ABS(il)««2-4*AA«(C-0,S«R3 
TYPL 4»DISC,AA,t»,C 
XK(ulSC.LT.O.w) GUTU 44 
ALPrtAl8(-b-SQRT(niSC))/(2TAA) 


»?-77 


ALPMA2sC-i)4>bOKT(OXSC) J/C2*AA) 

MU8 1 ♦ A liUG ( A L PH A ;< ) 1 1 . 0 
GUTU 33 

44 ACC6.pt*, NHJ 

33 ClJNTThni: 

nu ^5 tbi 

XA(I)s<».« 

OU 2f> iJsl,M 
Ar(l,J)aA(.1,13 

XX(ji3aXx(X) + A‘i'(i,JJ*D(J) 

2b CON'aTNUX: 

XX(X3sbXa(I)/ 81GAA 
25 CuNTlNne 

OU 27 Ts] 

XCXJau.O 

OU 2R Jsl,N 

SUHbu.O 

DU 29 Kb ],«1 

SUMbSUM+ATCI ,n)*A(K,J) 

29 CURTlAUe 

AT A X I , ) aSllA /SI 6H A 
X(1}BX(I)+AIA(I,J) 

2b CUNIlNflC 

27 COUTINUL 

ACCfcPT*, MOITjIACK 
DU JO iBlfN 
p(i)«6:xp(ixu*t.o) 

30 ruRxibiiis: 

KA«0 

LAMUbu.U 

77 CAllii CAbQF 

IF(ABS(uP}.liR«ePS3GUTU 99 
IFCVir.GT.O.OGQTn 31 
LAMusLAMO-InCR 
GuTO 32 

31 TiAPH«LAA0-l-IbC8 

32 CUNTTaUiL 
KKoKXfl 

IF(KR.GT.H01T)GUTU 999 
CATiLi SOLVe: 

DU 4C TslfK 
F(IJ»FF(IJ 
40 CONTI NUC 

LAMUbIjAmX 
GuTU 77 

99 CUNTInUL 

SUHsl'.O 
DU 50 TstfN 
SOMBSUMfFCI) 

PICiX}BiPiX(P(T)+0.5) 

50 CONTiMUfc: 

TIPS*, SUM, T 

HRTTGC45,*) (PICC X),Ib 1,N3,KK, LAMA, bAHO, SUM, T 
GOTO 55 





9«^9 TYPE lll,KK,QF,LAMQ,f.AMd 

111 PORdAl'(lX,*riO OF 1TERAT10M8 EXCERUEU*,Ib,iPi5.61 

C AFTe.P THIS NU. OF IT£RAT10«S AkR RAC£RUEu, PUE^iH VftLJR OF \itO 

C TACKCnEnT rOiCR ARld FRO UEPENDlNti ON TilE RATbl UP CO'^^RkG^WCT;. 

no 123 KBlfW 
piijBFPaj 
123 coNTiNne 

OOTO 77 
55 STOP 

B«0 

C ——————————— ™—— 

C SUPkQUTiNG CALQP CALCULATbS THE VALUE OP iPF-l/<t:4j .iPi'RK iilACM 

C 'iTRKATtUN Tu ACHIEVE OESImEU COrtVERiiEnCE. 

C 

SIIRKOUTINE CALOF 

nXMbNSIUN A(64,64j«D(b4)«F(64j,FPC64)rAF(b4j 
REAL LAHO,LAMN,Sl(iNA,UP#Su74,yx 
TrtTbGER N,M,L 

CUMrtQN/ARRAY8/A,D»P,FF,ATACb4,64),XXC64) 

CUMAOM/VAK/ttPfliAMUfLArtN 

COKAON/COASI/U r S f lif SitiSA 

QFBO.O 

DO 10 Isl,M 

SUP-0.0 

DC 20 Kal,S 

SUPs8UMfACl^K)«F(k) 

20 COHTIrlUE 

AFCl)aSUP-OCXl 

OF»«F+AFCI)vAFCI) 

10 coktinue 

0FB0,5*(CwF/SlGhA)-S) 
type 55b r of 

555 FURHAT(lX,'THe VALOE OP OF iS;*,P15,8j 

return 

END 

c 

C SUBKOuTiNb SOLVE IS MEANT FUR SuT.VInG A SxSlEn UP DlPFEHEI^TXAij 

C EgilATIONS BY OAUSS-SIEDEL IlEHATIVE SCHEPi!:. 

C 

SUHRUUTlNb SOLVb 

DXMERblUM Pl64»o4),BB(64)fFA(b4,64) 
real LAAn,r.AiN,STtiHA,OF,A,D,FrPF»ATArAX 
INTbGER M,M,L 

CClHHON/ARi<AYS/AC54r64) «D(t>4) »K(04) ,FP(64) ; ATA(64rt}4i ,XXl64) 
CUMHON/VAK/UFfLAHUfLAMN 
COPAON/COi«SX/ArS»L>SIGMA 
DO 10 Isl,N 
DO 10 Jsl,N 
FECX.J)sO.O 
10 CONTINUE 

DO 20 Tst,N 
• FKC1»1)«1«0/F(1) 

20 continue 

DU bO iBlfN 


60 


DO bO Jal,N 

P(I,JJsFKCl,J)-i>LAMO*AT4(IfJj 
conriHiiiu 

Do 70 £al,M 
SUMsO.O 
DO bO Jsl,4 
SUriaSUM+ATA( I , J) «K( J) 

80 CONTIMUe 

DbC i )sSUM* (2*LAM0 «LAMh )<ft .0 
70 CONTTHtlfc; 

DU OO Isl^N 

All C I ) bHB C I ) •«’ X X C 1 ) * ( L AMN-L AMQ ) 

90 CUK'VlNUi:: 

TeRMsu.U 

DU 110 ja>2fN 

TfiRh«TEKMtPCl,J)*r(J) 

110 CON‘J!lI«n& 

Ff (1 )s(bB(1)-TERM)/PC1,1) 
no 120 lai,A 

Tl!.RHa0,0 

DO l3o jal,l«l 

Tfe:RHal'ERM4-PCI»J)«PF(J) 

130 cuniinue 

PF(i)xRdCl)-TERH 
IF(Z.fi0«8)GUTO 199 
SUMad.O 
DO 140 

AUMaSUMtP(l,J}*F(J) 

1 40 CONTIHOE 

FF(1)*FF(I)-SUM 
199 PF(l)«FF(l)/PiI,I) 

120 continue 

RETURA 

RnD 

C 
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